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1  INTRODUCTION 


This  report  summarizes  the  research  performed  under  AFOSR  contract  49620-S9-C-0082, 
"Transform  Methods  for  Precision  i>'onlinear  Wave  Models  of  Flexible  Space  Structures."  The 
work  has  focused,  primarily,  on  exact  modelling  and  control  of  flexible  strucmres.  Two  modelling 
approaches  have  been  developed,  each  of  which  has  motivated  a  structural  control  methodology  in 
a  namral  way.  The  Transform  Element  Modelling  (TEM)  approach  uses  the  Laplace  transform  to 
obtain  exact,  frequency-domain  structural  element  models.  An  assembly  procedure  is  then  used  to 
create  an  exact  global  model.  Several  open-loop  control  algorithms  have  been  developed  based  on 
the  TEM  models  generated.  The  direct  PDE  approach  deals  specifically  with  the  underljing  time- 
domain  partial  differential  equation  describing  the  strucmral  behavior.  This  methodology  leads  to  a 
distributed  control  theory,  which  is  analogous  to  traditional  state-space  control.  These  iw’o  control 
algorithms  are  alternatives  to  conventional  structural  control  approaches.  They  are  designed  to 
alleviate  the  problems  associated  with  control-strucmre  interaction,  as  described  in  the  next  section. 

1.1  The  Control-Structure  Interaction  (CSI)  Problem 

The  requirements  for  many  military  and  civilian  structures  applications  both  in  space  and  on 
earth  call  for  the  use  of  large,  high  performance,  lightweight  strucmres.  In  most  cases,  the 
structural  weight  must  be  kept  as  small  as  possible  to  avoid  excessive  transport  costs.  However, 
the  flexibility  associated  with  large,  lightweight  structures  increases  the  likelihood  of  troublesome 
strucmral  behavior.  The  vibrational  modes  generally  begin  at  low  frequency,  and  are  excited  by 
disturbances.  Potential  sources  of  these  dismrbances  include  rotating  machinery  for  terrestrial 
applications  and  attimde  control,  antenna  retargeting  and  payload  shifting  for  space  applications. 
Often,  the  structure,  disturbance  and  control  bandwidths  arc  close  or  overlapping,  causing 
undesirable  vibration  to  propagate  through  the  structure.  Tnis  simation  is  described  in  Fig.  1-1.  In 
systems  with  stringent  pointing  requirements  (such  as  space-based  telescopes,  interferometers, 
lasers,  etc.),  this  can  severly  degrade  performance.  Similarly,  for  systems  requiring  pilot  isolation 
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Hg.  1-1;  Dcscripiioii  of  the  CSI  jiroblcin.  liiicracliun  occurs  when  the  siniclural  modes  lie 
witliin  the  bandwidth  of  the  control  system  and  disturbances. 


or  flutter  control,  performance  degradation  occurs  if  unwanted  vibration  is  present.  When  these 
problems  arise,  some  sort  of  vibration  suppression,  either  active  or  passive,  is  required. 

Passive  damping  techniques  alone  are  usually  insufficient  to  meet  performance 
requirements.  In  order  to  achieve  significant  modal  damping  (on  the  order  of  50%,  say),  an 
unacceptably  large  mass  of  damping  material  must  be  added  to  the  structure.  As  a  general  rule, 
passive  damping  can  only  provide  about  5%  to  10%  damping  for  a  5%  increase  in  strucmral  mass 
(See,  for  example,  Plunkett  (1970)).  Consequently,  these  techniques  provide  low  levels  of 
vibration  suppression  and  are  well  suited  for  addressing  steady-state  di.tturbances  arising  from  on¬ 
board  or  environmental  sources.  Active  control  techniques  are  required  for  suppressing  transient 
behavior  and  large,  dismrbance-induced  structural  responses. 

Large,  b'ghtweight  structures  have  three  basic  characteristics  which  make  active  control 
design  difficult.  First,  the  structures  are  difficult  to  model  precisely.  The  suucmral  dviiamics 
mathematical  models  obtained  from  finite  element  techniques  only  approximately  analyze  high 
frequency  behavior.  Typically,  only  the  first  few  modes  are  known  to  any  degree  of  accuracy.  As 
a  result,  the  active  control  design  must  be  extremely  robust.  Unfortunately,  robustness  often  leads 
to  V  ery  conserv  ative  designs  w  hich  sacrifice  performance  for  stability.  Second,  these  structures 
are  modally  dense  and  lightly  damped.  This  makes  the  closed-loop  system  extremely  sensitive  to 
para.me;er  v  ariarions,  and  often  leads  to  instability.  Third,  the  underlying  dynamics  are  of  infinite 
order  Asa  result,  traditional  full-order  linear  quadratic  regulator  (LQR)  and  linear  quadratic 
gaussian  (LQG)  control  designs  are  not  dtecily  applicable  to  these  types  of  systems. 

These  structural  charaaeristics  lead  to  a  phenomenon  described  by  Balas  (1978)  as 
"spillover,"  a.  d  can  be  explained  as  follows.  In  a  typical  control  design  procedure,  the 
tpp'cximatc  fimie  element  model  is  first  truncated  to  include  o.ily  those  modes  which  are  known  to 
a  good  degree  cf  accuncy  This  becomes  the  evaluation  model,  against  which  various  control 
dtiigta  art  j-dged  Typicady,  the  evaluation  model  is  of  too  high  an  order  to  achieve  an 
e-i.  -iCcrt  A*  a  miult,  a  reduced  order  controller  is  designed,  either  directly  from 

. «  e-fcrtrtv*.  ,  1985,), «  based  on  a  furfrier  truncation  of  the  evaluation  model 
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(Kosut  (1970),  Yousuff  (1984)).  In  eiihcr  case,  the  action  of  the  controller  excites  all  the  inodes  of 
the  evaluation  model  to  some  degree.  This  is  referred  to  as  control  spillover.  Likewise,  the 
sensors  associated  with  the  controller  also  respond  to  the  modes  truncated  from  the  evaluauon 
model,  leading  to  observation  spillover.  Additional  control  and  observation  spillover  occurs  when 
the  reduced  order  controller  is  applied  to  the  actual  infim’te  order  structure. 

1.2  The  Need  for  Better  Models  and  Control  Designs 

Chinrently,  the  methodology  generally  employed  in  most  structural  control  problems  is  the 
High-Authority  ControUer/Low-Auihoiity  Controller  (HA C/LAC)  approach,  described  in  Gupta 
(1981).  This  is  a  heirarchical  control  architecture  that  addresses  many  of  the  issues  presented  in 
the  presnous  section.  The  design  procedure  usually  involves  three  distinct  procedures.  In  addition 
to  developing  the  HAC  and  LAC  active  control  systems,  passive  damping  augmentation  is  usually 
designed.  In  many  cases,  prefiltering  of  command  inputs  is  also  required  to  minimize  excitation  of 
saucntral  modes.  The  modal-based  command  shaping  method  developed  by  Singer  (1990)  is  an 
example  of  the  prefiltering  concept. 

A  passive  damping  treatment  is  almost  always  required  in  the  control  design  of  infinite 
order,  lightly  damped  systems.  The  closed-loop  system  is  guaranteed  to  be  unstable  for  undamped 
infinite  order  systems  and  any  physically  realizable  controller  (i.e.,  any  controller  with  some 
amount  of  phase  lag).  This  phenomenon  is  explained  in  Fig.  1-2.  Furthermore,  even  slightly 
damped  systems  can  be  made  unstable  by  increasing  the  gain  of  the  controller  sufficiently.  As  a 
result,  passive  damping  is  usually  required  for  all  high  performance  structures.  For  some 
apph'cations,  adequate  passive  damping  may  be  inherently  present  in  the  struemre,  due  to  material 
friction,  viscoelasticity,  and/or  joint  hysteresis.  High  performance  applications  will  likely  require 
careful  tailoring  of  the  passive  design  in  order  to  meet  performance  goals.  By  shifting  the  open- 
loop  poles  of  the  system  into  the  left  half-plane,  passive  damping  has  the  added  oenefit  of 
desensitizing  the  controller  to  modelling  errors,  thereby  increasing  robustness.  Though  attractive 
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Fig.  1-2:  Collocated  rate  feedback  for  an  undamped  system.  The  dynamics  of  the  contro’ler, 
r‘presented  here  by  a  pole  on  the  real  axis,  alter  the  angles  of  departure  of  the  higher  frequency 
loci.  Instability  first  occurs  in  the  locus  with  approximate  radius  corresponding  to  the  radius  of  the 
first  order  pole  of  tlie  controller. 


Fig.  1-3:  Collocated  rate  feedback  for  a  damped  system.  Thepassivedamping,  inthiscasc,  is 
sufficient  to  prevent  instability,  and  the  damping  in  the  lower  modes  can  be  increased  significantly, 
In  an  actual  LAC  design,  many  actuators  and  sensors  are  used,  which  enables  greater  damping  of 
higher  modes. 
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for  practical  and  theoretical  reasons,  mass  penalo'es  place  an  upper  limit  to  the  amount  of  passive 
damping  that  can  be  implemented. 

The  LAC  (also  called  active  damping  augmentation)  is  usually  an  ad  hoe  design,  as  the 
primaiy  objecn’ve  is  to  achieve  robust  control  with  a  large  number  of  simple  controllers  (Aubrun 
(1980)).  Typically,  collocated  rate  feedback  is  employed.  The  effect  of  this  form  of  feedback  on 
the  modes  of  the  strucnnal  system  is  shown  in  Fig.  1-3.  The  main  objective  of  active  damping 
augmentation  is,  as  the  name  implies,  to  increase  modal  damping  significantly,  so  that  the  HAC 
does  not  destabilize  the  system  in  the  presence  of  modelling  errors. 

The  primaiy  design  objectives  are  accomplished  wa  the  HAC  It  is  usually  a  dj-namic 
compensator  of  high  order,  and  utilizes  information  from  sensors  located  throughout  the  structure. 
.Multiple  actuation  is  also  commonplace.  These  acnrators  and  sensors  need  not  be  distinct  from 
those  used  in  active  damping  augmentation.  As  the  name  implies,  the  HAC  exens  high  gain 
control  on  the  strucmre,  monng  the  closed-loop  pole  locations  considerably. 

The  current  state  of  the  an  in  HAOL.AC  design  suffers  from  several  serious  deficiencies. 
The  first  and  most  important  is  the  fidelity  of  the  underlying  structural  model.  The  traditional  firate 
element  modelling  approach  is  incapable  of  recovering  the  Wgh  frequency  •dj’namies  of  the 
structure,  unless  an  extremely  fine  discretization  is  utilized,  which  is  usually  computationally 
unacceptable.  This  limits  the  accuracy  of  the  evaluation  model  and  its  utility  in  compensator 
design.  Furthermore,  modelling  of  damping  mechanisms,  such  as  viscoelasticity,  is  difficult  to 
formulate  in  a  finite  element  environment.  The  damping  is  usually  assumed  to  be  model,  and  the 
actual  values  for  the  damping  ratios  are  determined  quite  arbitrarily. 

Clearly,  tlien,  a  more  accurate  structural  modelling  approach  would  be  extremely  beneficial. 
It  particular,  it  is  desirable  to  develop  a  modeUing  methodology  that  avoids  modal  truncation  and 
spatial  discretization  altogether.  This  is  the  basis  for  Partial  Differential  Equation  (PDE)  modelling 
approaches.  The  equations  describing  the  structure  are  kept  intact,  and  remain  iriathemarically 
exact  at  all  frequencies.  The  issue  of  modelling  error  is  then  reduced  to  knowledge  of  the  physical 
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parameters'  for  the  system  and  the  actual  choice  of  the  mathematical  representations  of  the  structural 
elements. 

Another  limitation  of  the  HAQLAC  methodology  is  in  the  active  damping  augmentation, 
which  is,  for  the  most  part,  ad  hoc.  Some  systematic  approaches  have  been  developed,  but  have 
limited  applicability.  For  example,  MacMartin  (1990)  describes  a  method  for  minimizing  the 
power  imparted  to  the  structure  at  an  interface  using  techniques.  However,  only  the  near  field 
effects  are  considered.  Refiections  of  the  dlsnirbances  at  other  boundaries  of  the  struemre  are  not 
modelled.  This  represents  a  worst  case  design,  with  the  assumption  that  nothing  is  known  about 
the  struemre  except  in  the  immediate  xicinity  of  the  controller.  This  approach  is  usually  too 
conservative  for  practical  applications.  Miller  (1991)  uses  a  wave  propagation  approach  to  design 
controllers  that  absorb  power  at  structural  interfaces.  Again,  the  lack  of  a  far  field  model  results  in 
a  consert-ative  design.  Furthermore,  neither  of  these  designs  is  applicable  to  applying  control  in 
the  interior  of  a  strucmral  element.  A  more  systematic  LAC  design  methodology  with  less 
conservatism  is  clearly  needed. 

Current  methods  of  designing  the  high-authority  controller  are  t>’pically  overconservative, 
as  the  closed-loop  systems  designed  from  the  truncated  model  tend  to  be  sensitive  to  parameter 
uncenainty.  Tjpically,  only  the  first  few  modes  are  used  Ln  the  design.  The  remaining  modes  are 
considered  urjeh'able,  and  contribute  to  modelling  error.  A  less  conservative  design  (and  therefore 
one  with  bener  nominal  performance)  would  arise  from  an  exact  strucmral  model,  as  the  modelling 
errors  are  smaller  and  are  more  easily  characterized. 

The  problems  discussed  above  suggest  that  new  modelling  approaches  may  lead  to  new 
control  design  methodologies  in  a  natural  way.  By  using  PDE  models,  mathematically  exact 
representations  of  structural  systems  are  available.  Control  designs  based  on  these  models  will 
account  for  all  modes  of  the  structure,  so  that  the  design  need  not  be  overly  conservative.  The 
achievement  of  this  design  methodology  would  significantly  enhance  the  overall  perfoimance  of 
large,  flexible  struemres  in  the  presence  of  disturbances. 


7 


1.3  Summary  of  First  Year  Research 

The  research  conducted  in  the  first  year  of  this  program  focused  on  developing  exact 
structure  models  that  preserve  the  vi’ave-like  eharacterisdcs  of  snuctural  disturbance  propagation. 
The  primary  emphasis  was  placed  on  developing  the  TEM  methodology  and  validating  it  with 
several  structural  evaluation  models.  Of  particular  concern  in  the  first  year  were  numerical  issues 
related  to  the  accuracy  of  the  model.  For  low  and  high  frequencies,  numerical  overflow  and 
roundoff  errors  occured  when  computing  structural  responses.  These  issues  were  resolved  by 
using  low-  and  high-frequency  asymptotic  approximations  where  required.  Also  developed  was 
an  open-loop,  optimal  control  algorithm  for  linear  and  small  angle  slews  of  flexible  structures.  The 
details  of  these  developments  can  be  found  in  Lupi  (1990),  but  are  also  presented  in  Chapters  2 
and  4  for  completeness.  The  application  of  the  TEM  approach  to  two-dimensional  structural 
elements  was  also  addressed.  In  Lupi  (1991b),  a  plane  stress  element  was  developed  that  utilized  a 
frequency-dependant  version  of  the  Airy  stress  function.  The  results  of  this  study  are  also 
presented  in  Section  2.4,2. 

1.4  Oterview  of  Report 

This  report  summarizes  developments  in  structural  modelling  and  control  design  from  the 
past  year  as  well  as  reviewing  the  previous  year  s  research.  Chapter  2  introduces  the  TEM 
methodology  and  its  applications  to  structural  analysis  problems.  The  extension  of  the  TEM 
approach  to  multibody  systems  is  addressed  in  Chapter  3.  Open-loop  conuol  designs,  based  on 
TEM  models,  are  presented  in  Chapter  4.  The  extension  to  closed-loop  control  is  also  discussed. 
Chapter  5  develops  the  direct  PDE  modelling  approach  for  one-  and  two-dimensional  strucniral 
elements.  The  direct  simulation  of  a  simple  one-dimensional  system  is  also  presented.  The  direct 
PDE  modelling  approach  leads  naturally  to  a  distributed  control  theory,  which  is  developed  in 
Chapter  6.  The  possibility  of  a  hybrid  modelling  approach  for  HA  C/LAC  control  designs  is 
discussed  in  Chapter  7.  Finally,  conclusions  and  recommendations  are  presented  in  Chapter  8. 


2  THE  TEM  MODELLING  APPROACH 


In  this  chapter,  we  introduce  the  Transform  Element  Modelling  (TEM)  method  of  strucniral 
analysis.  This  approach  begins  with  the  partial  differential  equations  describmg  the  dynamics  of 
the  inditidual  elements  that  describe  the  strucntre.  The  Laplace  transform  is  then  utiliaed  to 
express  these  equations  in  the  frequency-domain.  Exact,  frequency-dependant  stiffness  matrices, 
which  relate  a  set  of  generalized  forces  to  a  set  of  generalized  displacements  at  the  element 
boundaries,  are  then  derived.  These  stiffness  matrices  are  then  assembled  to  form  a  global  model, 
in  a  manner  similar  to  the  traditional  finite  element  assembly  procedure.  The  advantages  of  this 
approach  are  described  in  Section  2.1.  The  procedure  for  converting  frequency-domain  data  back 
into  the  time-domain  is  described  in  Section  2.2.  The  general  theory  for  one-dimensional 
elements,  and  some  simple  representative  examples,  are  presented  in  Section  2.3,  while  Section 
2.4  applies  the  TEM  methodology  to  two-dimensional  structural  element  models.  The  assembly 
procedure  is  desenbed  in  Section  2.5,  and  several  applications  of  the  TEM  approach  are  presented 
in  Section  2.6. 

2.1  Advantages  of  the  TEM  Method 

The  TEM  approach  has  several  important  advantages  over  the  traditional  Finite  Element 
Method  (FEM).  as  measured  in  temis  of  numerical  accuracy  and  computational  efficiency.  The 
dynamics  of  the  one-dimensional  elements  that  comprise  the  structure  to  be  analyzed  are 
represented  exactly  in  the  TEM  approach.  This  is  made  possible  by  the  Laplace  transform 
operation,  which  converts  the  time-domain  partial  differential  equation  of  a  structural  element  into  a 
frequency-dependant  ordinary  differential  equation.  As  a  result,  analytical  general  solutions  to 
these  equations  are  available  for  most  common  element  models,  such  as  Bernoulli- Euler  and 
Timoshenko  beams,  and  axial  and  torsional  rods.  More  complicated  elements  can  be  handled  by 
special  numerical  modelling  approaches.  This  is  in  marked  contrast  to  FEM  modelling,  where  each 
element  is  divided  into  smaller  subsections,  each  of  which  is  constrained  to  deform  with  a  finite 
number  of  degrees  of  freedom.  The  interpolation  functions  associated  with  the  deformational 


degrees  of  freedom  satisfy  the  underlying  differential  equation  exactly  for  the  stanc  case,  but  only 
approximate  the  exact  solution  the  dynamic  case.  Consequently,  for  dynamics  problems,  the  FEM 
analysis  can  only  yield  approximate  results.  Conversely,  the  TEM  approach,  which  utilizes 
frequency-dependent  (generally  transcendental)  interpolation  functions,  is  mathematically  exaa  at 
all  frequencies.  Furthermore,  since  only  one  mathematical  element  is  required  for  each  physical 
strucnual  element,  the  TEM  approach  is  far  superior  to  the  FEM  methodology  in  terms  of 
computational  speed. 

For  two-dimensional  strucmral  elements,  general  exact  bJ'  tions  are  not  available,  but  the 
TEM  methodology  makes  it  possible  to  approximate  the  solutions  in  terms  of  finite  series  of 
displacement  functions.  Each  of  these  functions  satisfies  the  underlying  differential  equation  in  the 
interior  of  the  element  exactly,  and  approximations  are  made  only  at  the  boundaries.  A  comparison 
between  the  TEM  and  FEM  methods  in  terms  of  speed  and  accuracy  for  two-dimensional  elements 
has  not  yet  been  attempted. 

Another  imponant  advantage  of  the  TEM  approach  is  its  ability  to  incorporate  general  linear 
viscoelastic  damping  models  in  a  straightforward  manner.  Using  the  correspondance  principle,  as 
described  by  Hughes  (1 989),  the  physical  parameter  of  interest  is  simply  replaced  with  a 
frequency-dependent  counterpart.  For  example,  the  Voigt  damping  mechanism  is  easily  expressed 
by 

E(s)  =  [l+<iv^]Eo  (2.1) 


where  E  is  the  modulus  of  elasticity  of  the  material.  Eg  is  its  static  value,  s  is  the  (generally 
complex)  frequency,  cOy  is  a  characteristic  frequency,  and  cty  is  an  empirically  determined 
nondimensional  parameter.  A  more  general  damping  model,  suggested  by  Hughes,  is 


E(s) 


Oi 


s2+2(;,W,5+(0j 


Eo 


(2.2) 
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where,  for  each  value  of q  is  a  characterisnc  scaling  factor,  cO;  is  a  characteristic  frequency,  and 
C;  is  a  c’'aracteristic  damping  ratio,  all  of  which  are  empirically  determined.  However,  the 
advantage  of  the  TEM  approach  lies  in  its  ability  to  model  damping  mechanisms  of  infinite  order, 
such  as  the  fracdonal  derivative  models  discussed  by  Bagley  (1983).  Such  models  require  the  use 
of  fractional  calculus  techniques  when  employed  hr  the  time  domian,  but  are  easily  cast  in  the 
frequency-domain  as  fractional  powers  of  the  complex  frequency.  For  example,  the  square  root 
damping  model  is  wrinen  as 


where  Oj  is  a  characteristic  frequency  and  q  is  empirically  determined.  Thus,  since  all  linear 
s  iscoelasric  damping  models  have  frequency-domain  representations,  any  model  can  be  used  in  the 
TEM  fonnulation. 

Yet  another  advantage  of  the  TEM  approach  is  the  ability  to  take  derivatives  of  time-domaut 
functions.  All  that  is  required  is  multipUcadon  of  the  function  by  the  complex  frequency  r’ariable. 
Similarly,  additional  multiplications  by  s  jield  higher  order  derivatives.  Thus,  given  a  set  of 
frequency-domain  dau  representing  a  time-domain  response,  it  is  a  simple'matter  of  multiplying 
the  data  by  the  complex  frequency  before  invoking  the  inverse  Laplace  transfoim  to  obtain  the 
derivative  of  the  response. 

2.2  Inverse  Laplace  transform  algorithm 

In  any  frequency-domain  modelling  approach,  it  is  of  paramount  importance  to  have  the 
ability  to  convert  data  back  into  the  time-domain  in  a  computationally  efficient  manner.  This  is  the 
basis  for  many  inverse  Laplace  transform  algorithms,  such  as  the  method  of  expansion  by 
Laguerre  functions  as  described  by  Weeks  (1966)  and  Wing  (1967).  However,  the  most 
straightforward,  stable  and  accurate  method  for  general  functions  appears  to  be  the  numerical 
approach  of  VVilcox  (1978).  It  is  this  method  that  has  been  used  exclusively  in  this  research,  and  it 
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therefore  deserves  tnennon  here.  A  detailed  comparison  of  several  other  approaches  may  be  found 
in  Davies  (1979). 

The  Laplace  transfonn  pair  is  expressed  as 


i(s)  =  Jg(t)e“dt 
g(t)  =  Jg-(s)e‘'ds 


(2.4a, b) 


where  C  is  a  closed  contour  which  encloses  all  the  singularities  of  g(s).  If  we  assume  that  g(s) 
has  no  poles  with  real  parts  greater  than  ot,  where  a  is  a  small  positive  number,  then  the  integration 
contour  can  be  taken  as  the  Bromwich  contour,  which  is  shown  in  Fig.  2- 1 .  Furthermore, 
assuming  that  the  semicircle  pan  of  the  contour  does  not  contribute  to  the  integral,  the  inverse 
transform  reduces  to 

oo 

g(t)  =  ^  |g'(s)e'“'do),  s  =  a+jco  (2.5) 


Since  the  path  of  integration  is  displaced  to  the  right  of  the  imaginary  axis,  marginally  stable  and 
slightly  unstable  functions  can  be  invened.  Because  g(t)  is  assumed  to  be  real-valued,  it  follows 
from  Equation  (2.4a)  that  g(s*)  =  g(s)*.  Asaresult,  Equation  (2.5)  reduces 'o 

c« 

g(t)  =  e"  Re  {  J  J  g' (s)  ej“'  do) }  (2.6) 

The  numerical  computation  of  g(i)  involves  calculating  g(s)  at  N  evenly-spaced  complex 
frequencies  along  the  Bromwich  coi.tour,  as  shown  in  Fig.  2-2.  These  values  are  given  by 

U)c  =  (2k-i-l)Ao),  k  =  0 . N-1  (2.7) 

The  algorithm  yields  2N  values  of  g(t)  at  evenly-spaced  time  intervals,  given  by 

tm=^T.  m  =  0 . 2N-i  (2.8) 
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where  the  time  interval  extends  from  t=0  to  t=T.  Wilcox  (1978)  shows  th ut,  for  redprocity,  the 
following  relationship  must  hold: 


Am  =  I 

Hence,  Equation  (2.6)  can  be  approximated  by  the  midpoint  rule  as  follows: 

N-1 

k=() 

Factoring  out  constant  terms  from  within  the  summation  and  simplifying  yields 

N-I 

S(tJ  =  e“'--"Re  {  Z'^”} 

k=0 

where 

Z  =  ei^--'-’^ 

If  w’e  now  define  gj.  and  g^,  by 

gk  =  i(Sk)Z(=^-»/^ ,  g(t„,)  =  e“™Re  {  ^Z’^'^g^  } 

Then  Eq.  (2.1 1)  becomes 

N-I 

Sm  =  ZskZ’-’" 
k=0 


(2.9) 


(2.10) 


(2.11) 


(2.12) 


(2.!;5a,b) 


(2.14) 


For  computarional  efficiency,  it  is  useful  to  write  Eq.  (2.14)  in  a  form  amenable  to  fast  Fourier 
transform  techniques.  This  is  accomplished  by  separating  the  time-domain  samples  into  even  and 
odd  sets.  Thus, 

!,"  “  I  n=0,...,N-l  (2.15) 

?n  =  g2n+l  J 
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and  the  inverse  transfora  is  applied  nvice,  yielding 
N-l 

in  =  EgkW^. 

lc=0 


N-l 

gn  =  SgicW‘=" 

k=0 


where 


gk  < 


w  = 


(2.1 6a, b) 


(2.17a, b) 


Thus,  given  a  series  of  samples,  i(sO,  in  she  frequency-domain,  the  algorithm  is  as 
follows:  Use  Eq.  2.13a  lo  obtain  gj;  ar,d  Eq.  2.17a  to  obtain  gj;.  Next,  apply  the  fast  Fourier 
transform,  as  described  by  Cooley  (1970),  to  obtain  and  gn.  Upon  reordering  the  data,  which 
yields  g^,  use  Eq.  2.13b  to  obtain  g(tm).  Wilcox  (1978)  shows  that,  as  a  general  rule  of  thumb,  it 
is  best  to  take  a=2rJT.  Using  this  rule,  the  author  has  determined  the  algorithm  to  be  accurate  to 
0.1%  for  the  first  75%  of  the  simulation  for  all  test  functions,  with  some  deterioration  occuring 
after  this  time.  This  can  be  o\’ercome  by  increasing  the  simulation  time  slightly  and  discarding  the 
later  data 

In  eases  where  the  time-domain  response  contains  step  discontinuities,  it  is  useful  to  scale 
the  frequency-domain  data  by  a  Gibbs’  oscillation  suppression  factor,  given  by 

sin[(2k-fl);t/2N] 


fk  = 


(2k+l)rt/2N 


(2.18) 


This  has  the  equivalent  effect  in  the  time-domain  of  passing  the  signal  through  a  fmite-time 
integrator  with  time  constant  equal  to  the  time  between  samples,  as  explained  by  Lancaos  (1957). 
Consequently,  this  scaling  does  not  affect  the  response  where  it  is  continuous  in  time,  while  it 
reduces  Gibbs'  oscillations  considerably  at  discontinuities  (at  the  expense  of  slightly  increased  rise 
time).  Fig.  2-3  shows  some  time-domain  responses  generated  by  the  inverse  Laplace  transform 
algorithm,  both  with  and  without  the  Gibbs  suppression  faett  .  The  favorable  effect  of  the  scaling 
is  obrious.  Thus,  a  numerically  robust  inversion  algorithm  has  been  presented.  This  algorithm 
has  been  used  extensively  in  this  reasearch. 
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2.3  0ne-dimensicn2l  element  models 

We  now  proceed  to  develop  the  TEM  formulanon  for  one-dimensional  elements.  The 
general  equation  of  isoncn  for  a  one-dimensional  structural  clement  can  be  wrinen  as 

ku4[v()t.0]+l:T''(x.')  =  Wx.O.  xe  [O.L] ,  ts  [0,~)  (2.19) 

»he«  v(x,t)  IS  a  generahaed  displacement,  fjj(x,t)  is  a  generaliaed  dismbuted  forcing  function,  L, 
is  a  linear  spadal  differential  operator  of  order  n,  and  (')  denotes  differennatioa  with  respect  to 
dme.  The  constants  ky  and  kj  are  physical  parameters  related  to  the  internal  potential  and  Idneric 
energy  of  defo.mation  of  the  strucniral  element,  respectively.  They  can  be  thought  of  as 
generalized  stiffness  and  mass  parameters.  The  boundary  conditions  are  as  yet  unspecified.  The 
equation  of  morion  is  such  that  all  ri,le\  ant  internal  states  (force  resultants,  moments,  etc.)  can  be 
obtained  via  spatial  differential  operations  on  v(x,t). 

2.3.1  Stiffness  matrix  formulation 

For  a  one-dimensional  element,  it  is  possible  to  obtain  an  exact,  frequency-dependent 
stiffness  matnx  relating  generaliaed  boundary  forces  and  displacements.  This  is  possible  because 
tiie  Laplace  transform  operation  converts  the  partial  differentia]  equation  into  an  ordinary 
differential  equation,  whose  solution  can  be  expressed  analytically  (usually  in  the  form  of 
transcendental  functions  of  die  complex  frequency).  Taking  the  Laplace  transform  of  Eq.  (2  19) 
leads  to 

I-,[v(x.s)] -f^s^vfx.s)  =  j^[fd(x.s) +  kxVo(x)  +k7SVo(x)]  (2.20) 

where  s  is  the  (generally  complex-calued)  Laplace  variable,  (')  denotes  the  transform  of  a  function, 
and  vo(x)  and  vq(,x)  represent  the  initial  conditions.  From  this  point  on,  the  overbar  on 
transformed  functions  will  be  assumed,  so  as  to  simplify  the  notation.  Also,  Jie  right  hand  side  of 
Eq.  (2.20)  will  be  lumped  into  a  single  function,  fd(x,s),  in  the  frequency-domain.  This  leads  to 
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(2.21) 


la['’(x.s)]+^s2v(x,s)  =  fd(x,s) 

Ws  wish  to  express  the  solution  to  the  preceeding  eqnarion  in  terms  of  generalized 
displacements  and  forces  at  the  boundaries  of  the  structurul  element  This  facilitates  the  assembly 
of  these  strucniral  elements  into  a  general  structural  model,  as  will  be  discussed  later.  Thus,  the 
genera!  soluuon  is  expressed  as 

L 

v(x.s)  =  v„(x,s)'’'a(s)  +  J'vp(x,^,s)  f(i(ts)d5  (2.22) 

where  vj^fx.s)  is  an  n-vecior  of  homogeneous  solutions  to  Eq.  (2.2 1 )  and  a(s)  is  an  n-vector  of 
arbitrary  constants.  The  Green's  function,  Vp(x,£,s),  comesponds  to  the  operator  +  (IttA'i))  s2 
and  satisfies  the  essentia!  homogeneous  boundary  conditions  at  both  boundaries. 

We  must  now  express  the  generalized  boundary  displacements,  vr(s),  and  generalized 
boundary  forces,  q(s),  in  terms  of  a(s).  Since  knowledge  of  v(x,s)  implies  knowledge  of  the 
entire  state  throughout  the  element,  these  boundary  conditions  are  ob’,ained  by  simply  evaluating 
v(x,s)  (and  its  derivatives)  as  given  in  Eq.  (2.22)  at  the  boundaries.  This  leads  to  linear 
relationships  between  the  boundary  states  and  the  arbitrary  constants,  which  can  be  expressed  as 

w(s)  =  'i'(s)a(s)  (2.23) 

and 

q(s)  = '?(s)a(s)  +  qo(s)  (2.24) 

where  '^(s)  and  y(s)  are  n-by-n  matrices,  and  q(j(s)  is  an  n-vector  arising  from  the  integral  term  on 
the  right  hand  side  of  Eq.  (2.22).  Because  satisfies  Vp(x,£,s)  satisfies  (by  construction)  the 
homogeneous  essential  boundary  conditions,  there  is  no  term  in  Eq.  (2.23)  corresponding  to  q(j(s). 
Combining  Eqs.  (2.23)  and  (2.24)  leads  to  the  desired  relationship  between  the  boundary  forces 
and  displacements,  giver,  by 
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q(s)  =  K(s)w(s)  +  q(i(s) 


(2.25) 


where  K(s)  is  refentd  to  as  the  dynamic  stiffness  matrix,  and  is  given  by 

K(s)  =  ’?(s)['I'(s)]‘‘  (2.26) 

It  should  be  noted  that,  at  any  complex  frequency,  Eq.  (2.25)  is  mathemadcaUy  exact.  This 
is  in  contrast  to  traditional  finite  element  stiffness  matrices,  which  are  usually  derived  from  an 
approximate  solution  to  the  equation  of  morion  describing  the  dynamics  of  the  strucmral  element 

2.3.2  interpolation 

The  exact  representarion  of  the  structural  element  is  not  restricted  to  the  boundary  forces 
and  displacements.  In  addition,  the  internal  states  of  the  element  at  an  arbitrary  location  can  be 
computed  exactly.  This  is  accomplished  by  first  expressing  the  n-dimensional  internal  stale  vector, 
u(x,s),  in  terms  of  v(x,s): 

lu(x,s)  =  L,(v(x,$)]  (2.27) 

Here,  Lx  “s  an  n-dimensional  spatial  differential  operator  vector,  and  the  superscript  (■)  indicates 
that  the  elements  of  u(x,s)  are  e>  pressed  with  respect  to  an  inertial  frame.  -Making  use  of  Eq. 
(2.22),  we  obtain 

’u(x.s)  =  0(x,s)a(s)-i-Up(x,s)  (2.2S) 

where  the  following  definitions  have  been  employed: 

L 

<I>(x,s)  =  L,[vh(x,s)'’']  ,  Up{x,s)  =  j Lx[vp(x,^.s)]  fd(£.s)d^  (2.29) 

Using  Eq.  (2.23)  to  eliminate  a(s)  jields 

'u(x,s)  =  C)(x.s)['F(s)]'‘w(s)  +  Up(x,s)  (2.30) 
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Thus,  ihe  iniemal  stales  are  easily  expressed  in  terms  of  the  boundary  displacements. 

It  is  usually  more  desirable  to  express  the  internal  states  in  a  frame  fixed  to  one  of  the 
boundaries  of  the  strucnrral  element  This  is  useful,  for  example,  if  the  structure  is  undergoing  a 
rigid  monoa  Expressing  the  internal  states  wiili  reqrect  to  a  fiame  fixed  to  an  element  boundary 
would  then  indicate  the  amount  of  internal  strucnrral  deformation  only.  This  change  of  reference 
affects  the  internal  generalized  displacements  only,  since  the  internal  generalized  forces  are 
automatically  zero  for  rigid  morion.  The  general  linearized  relationship  between  the  state  vector 
expressed  in  the  wo  frames  is  given  by 

u(x,s)  =  fu{x,s)  -  0(x)  w(s)  (2-31) 


where  u(x,s)  is  the  state  vector  expressed  in  the  roonng  frame  and  ©(x)  is  an  n-by-n  matrix  which 
is  independant  of  frequency.  Collecting  the  prewous  two  equations  yields 

u(x,s)  =  r(x,s)w(s)  +  Up(x,s)  (2.32) 


where 


r(x,s) 


0(x,s)['i'(s)]'’-e(x)- 


(2.33) 


A  particular  element  of  the  internal  state  vector  is  then  given  by 

Ui(x.s)  =  vi(x,s)'’‘w(s)  +  Up|(x.s)  (2.34) 


Thus,  once  the  generalized  displacements  at  the  boundaries  are  known,  it  is  a  simple  maner  to 
obtain  the  internal  states.  Once  again,  the  formulation  is  exact,  and  no  modal  truncation  or  finite 
element  approximations  have  been  made. 
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2.3.3  Intema]  energy  fonnuJaoon 

A  useful  scalar  measure  of  the  state  of  defonnatfon  of  a  flexible  element  is  the  total  energy 
due  to  deformation.  It  is  particularly  useful  in  control  appUcadons,  as  it  is  a  quadratic  fiinaion  of 
defoimation  amplitude  and  is  therefore  well  suited  for  linear  quadratic  regulator  problems.  The 
internal  energy  within  a  particular  strucmral  element  is  obtained  by  integrating  the  sum  of  the 
potential  and  Idnenc  energy  densities  over  the  length  of  the  element.  This  leads  to  an  expression  of 
the  form 

L 

E(t)  =  5  J  [  kyuSfx.t)  +  k'Tblfx.t)]  d  X  (2.35) 

where  E(i)  is  the  internal  energy  at  time  t,  and  Uu(x,t)  and  u-j-fx.t)  are  components  of  the  internal 
state  vector  related  to  the  potential  and  kinetic  energies,  respectively.  Note  that  these  components 
are  now  expressed  in  the  time-domain.  In  order  that  the  internal  energy  be  independant  of  rigid 
morion,  it  is  imperative  that  UpCx,!)  be  expressed  with  respect  to  one  of  the  boundaries  of  the 
element,  as  described  in  the  prexious  section. 

The  terms  Uu(x,i)  and  u^fx.t)  are  new  e.xpressed  as  the  inverse  Laplace  transforms  of  tite 
comesponding  frequency-domain  functions.  This  is  accomplished  via  Eq.  (2.5)  and  leads  to 

L  «>  2  «  2 

E(0  =  (  J Uu(’t.s)«^“'tiw)  +  kT(  Js  UT(x.s)e^“'dco)  )  dx  (2.36) 

where 

k(0  =  (2-37) 

We  seek  an  exact  expression  for  the  total  energy  of  deformation.  It  is  therefore  necessary  to 
replace  the  squared  integrals  with  double  integrals,  so  that  the  order  of  integration  with  respea  to 
space  and  frequency  may  be  reversed.  This  makes  it  possible  to  perform  the  spatial  integration 
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analytically.  Thus,  by  wilting  each  inverse  transform  integral  next  to  itself,  using  different  dummy 
variables  of  integration,  and  grouping  the  integrals  together,  we  obtain 


L 

E(t)  =  Juu(X.Si)Uu(x.S2)ei<“>*"2)'d(D,dCD2 


+  hi"  J  j  s  1  S2  ux(x,si)UT{x,S2)e'^'*“^^'dcoj  dco2 }  dx  (2.38) 


where 

Sj  =  a+jtOj  (2.39) 

For  simplicity,  we  have  assumed  that  the  initial  conditions  are  zero  and  that  no  distributed  forcing 
occurs  in  the  interior  of  the  element.  Interchanging  the  order  of  the  spatial  and  frequency 
integrations  in  Eq.  (2.38)  yields 


E(t) 


^  0 


Si)uu(x.S2)dx 


+  kjS] S2J ux(x,si) uj(x,S2)dx  }  e'^“>*^^dcoj  dcoi  (2.'40) 

We  are  now  able  to  express  the  spatial  integrals  in  terms  of  the  boundary  displacements.  Making 
use  of  Eq.  (2.34)  leads  to 

L 

Jui(x,Si)Ui(x,S2)dx  =  w(si)'''’Ei(sj.S2)w(s2)  (2.41) 

where 
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(2.42) 


L 

Hi(si.S2)  =  jYi(x.Si)Y,(x,S2)''‘dx 

Since  Yu(x*sj)  and  Yr(x.Si)  are  expressed  analyncally,  the  matrices  Hu(si,S2)  and  Sj<S],S2)  can 
also  be  computed  exactly  at  each  frequency.  Finally,  substintting  Eq.  (2.41)  into  Eq.  (2.40)  )ie!ds 


E(t)  = 


Wi) 

2 


J  I w(si)'''E(sj,S2)  w(S2)eJ^“>'^“25'do)i  dt02 


(2.43) 


where 


— (S]iS2)  ~  l;tjEu(S],S2)  +  kxS)S2=T(S),S2)  (2.44) 


Thus,  given  the  boundary  displacements  in  the  frequency-domain,  the  energy  of 
deformation  is  computable  via  a  double  integral.  For  general  motions,  an  analytical  solution  does 
not  exist,  and  E(i)  must  be  co.mputed  numerically.  The  compulation  time  is  reduced  by  a  factor  of 
two  by  exploiting  the  following  symmetry  properdes  of  E: 


E(S2,Si)  =  z(Si.S2f 
E(SpS2)  =  Efsi.sj)’, 


(2.45b,c) 


UsLng  these  properties  makes  it  possible  to  reduce  the  integral  to 


+02)1 


E(t)  =  k(t)  J J  [w(si)'''S(si,S2)w(S2)e)^®>'"“ 

+  w(s;)^  Efsj.sj)  w(S2)*c)*“‘’“^]  dCOj  d(02 


(2.46) 


In  the  actual  implementation  of  tiiis  formula,  the  integrals  are  replaced  with  summations,  the  upper 
limits  are  replaced  with  finite  frequencies,  and  a  simple  midpoint  rule  algorithm  is  invoked. 
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2.3.4  Axial  rod  example 

We  first  consider  a  unifonn  rod  constrdned  to  deform  axially,  as  shown  in  Fig.  2-4.  For 
this  element,  the  simplest  model  that  describes  its  djuamic  behavior  is  the  wave  equation,  given  by 

-EA^v(x,t)  +  pAv(x,t)  =  fd(x,t)  (2.47) 


where  v  represents  the  axial  deflection  of  the  cross  section,  A  is  the  cross  sectional  are,  and  p  and 
E  are  the  material  density  and  modulus  of  elasticity,  respectively.  Implicit  in  this  model  is  the 
assumption  that  the  deformation  of  the  element  is  uniform  across  the  cross  section.  Also, 
Poisson's  ratio  effects  are  ignored.  The  internal  state  at  any  location,  x,  is  therefore  completely 
chaiacteriaed  by  two  components: 


u(x,s)  = 


(2.48) 


Here,  F  represents  the  net  force  resultant  within  the  rod.  The  generalized  boundary  forces  and 
displacements  are  just  these  quantities  evaluated  at  x=0  and  xssL: 


■-Ea|;v(0,5)' 

EA^v(!.,s). 


(2.49a, b) 


The  homogeneous  solution  vector  is  simply 

vk(x.s)^  =  P  =  A/f’s  (2.50) 


and  the  Green's  function  kernel  is 


Vp(x,E,s) 


X  S  ^ 

x  >  4 


(2.51) 


Equations  (2.49a,b)  and  (2.50)  can  be  combined  to  determine  the  dynamics  stiffness  matrix,  which 
is  given  by 
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Fig.  2-5:  Bemoulli-Euler  beam  model.  Planar  cross  sccnons  remain  planar  and  perpendicular  to 
defonned  beam  axis. 
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VM  _  EAB  fcosh  PL  -1  ■] 

“  sinbpL  L  -1  cosh  PL  J 


(2.52) 


As  expected,  the  elements  of  K  are  transcendental  functions  of  the  complex  frequency.  Also,  this 
stif&ess  matrix  reduces  to  the  finite  element  static  stiffness  matrix  for  the  axial  rod  in  the  limit  as  s 
approaches  aero.  The  effects  of  initial  conditions  and  distributed  forcing  are  computed  via 

L 

“o®  =  -life j  fd(x.s)dx  (2.53) 


where  the  interpolation  matrix  is  given  by 

r-,  ,  1  r  sinhP(L-x)  sinh  px"!  Pi  Ol 

“  sinh  PL  L- P  cosh  P(L-x)  pcoshpxj  "  Lo  Oj 


(2.54) 


The  iniemal  energy  matrices  are  quite  complex,  but  they  arc  nonetheless  expressible  in  an 
analytical  form.  The  kinetic  energy  kernel  can  be  expressed  as  the  sum  of  the  following  four 
matrices: 


Ht(Si.S2)  =  i^A(j3i)T 


^{e-(PfP2)L  JJ  ^[e.(pK32)L  ,J 


i-Mhf 


o' 

Lg[e-?i"-3)  0_ 


1  f  [e-P^‘'-l)' 

0  0. 


A(P2)  +  [o  o] 


A(P2) 


(2.55) 


while  the  potential  energy  kernel  is  given  by 
=t;(si.S2)  =  ||^A(Pi)’‘ 


^[e(Pl*?2)Lj3  ^[^(Pl.PzlL.ij 

prP2  KJ-P2 


A(p2) 


(2.56) 


where  the  following  definitions  have  been  made: 
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(2.57a) 


f2.57b) 
(2.57c) 

It  IS  clear  that,  for  e\  en  the  simplest  of  structural  element  models,  the  internal  energy  expressions 
become  quite  complex.  However,  this  is  to  be  expected,  as  the  energy  is  a  quadratic  function  of 
deflection  ampliwde,  and  represents  an  integrated  effect  over  the  entire  domain  of  the  element. 

The  preceeding  expressions  for  the  stiffness  and  interpolation  matrices  are  well  suited  for 
numerical  computation,  pro\  ided  that  the  complex  ffequen-y  at  which  they  are  es'aluated  is  neither 
too  large  nor  too  small  in  magnitude.  For  these  extreme  situadons,  numerical  accuracy  and 
o\erf]ow  errors  become  issues.  These  proble.ms  are  readily  handled  by  using  asymptotic 
approximadons  to  die  sdffness  and  interpolation  matrices.  The  approximations  are  obtained  by 
replacing  the  trigonometric  and  hyperbolic  functions  with  appropriate  series  expansions,  and 
truncating  higher  order  terms. 


2  3.5  Bemoulli-Euler  beam  example 

For  bending  elements,  the  simplest  model  is  die  Bemoulli-Euler  beam,  shown 
schematically  in  Fig.  2-5.  The  basic  assumptions  of  the  model  are  that  planar  cross  sections  of  the 
bea.m  remain  planar  and  normal  to  the  center-line  after  deformation,  and  that  differential  cross 
sections  ha\  e  negligible  rotary  inenia.  Under  these  assumptions,  the  equation  of  motion  becomes 

EI^  v(x.t)  -s  pA’v(x.t)  =  fd(x,t)  (2.5S) 

w  here  v  is  the  transverse  deflection,  E  and  p  are  the  modulus  of  elasticity  and  density  of  the 
material,  respectively,  and  A  and  I  are  the  aoss  sectional  area  and  moment  of  inenia,  respectively. 
TaJting  the  Laplace  transform,  w'c  obtain 
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|^v(x.s)-aVx.s)  =  fd(x,s).  a-*  =  .e|£ 


(2.59a.b) 


For  this  eltiaent,  the  structural  state  vector  has  four  eJements,  given  by 


~  1 

^  A 

a* 

=  MM  =  V(X.S) 

Ls(x.s)J  "3 


where  0  is  the  rotation  of  the  cross  section,  and  M  and  S  are  the  interna]  moment  and  shear 
resultants,  respectively.  The  generalized  boundary  displacements  and  forces  are  then 


,  r 

v(0,s)  ^  , 

_  6(0.0  _ 

-  s(L.i)  -  v(L,s) 


-s(o,sn 

~  S(L.s)  - 


EI^v(0,s) 

■El|^v(0,s) 


.El^^(L,s) 

EI^v(L,s)_ 


(2.61a,b) 


The  homogeneous  solua'on  x'cctor  contains  four  elements  as  well,  and  is  gi\’en  by 
vij(x,s)^  =  [e®*  e‘®*  eJ“*  e'j®’'] 


and  the  Green's  function  kernel  follow's  as: 


Vp(x,C,S)  = 


gl(c.s)  S'(ax)  g;(^.s)  c'(c>)  . 

(1  •  ch  ct)  ^  ^ 

£,(i.s)  S-(a(L.x))  *  g^a.f)C'(c(L-.)) 

4a6  (1  -  cb  ct)  >  9 


gj(c,s)  «  -  C*(a;)  +  ch  cos(a(L-5))  •  sh  sa)(a{L-5))  +  ct  cosb(a(L-5))  +  si  sinh(a.1.-5)) 

S2(i.f)  =  S\ai)  +  ch  sin(a(L-c))  •  sh  cos{a(L-4))  +  ct  sinh(a(L-4))  •  st  cosh(a(l.-i)) 
53(4, s)  =  •  C*(o(L-4))  +  ch  co5(o4)  •  sh  si!i(a4)  +  ct  cosh(ci4)  +  si  $inh(a4) 

S',(C.s)  =  S*(o(I--4))  +  ch  5in{a4)  ■  sh  cos(a4)  +  ci  sinh(o4)  •  st  cosh(ci4) 


and  the  following  trigonometric  definitions  have  been  made: 
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C'^(o5)  =  cosh(a^)  +  cos(a£)  C'(al)  =  cosh(a4)  -  cos(a^) 
S^Ca^)  =  sinh(a5)  +  sin(a4)  S'(o£)  =  sinh(o^)  -  sin(a4) 
ch  =  cosh(oL)  ct  =  cos(aL) 

sh  =  sinh(aL)  st  =  sin(aL) 


For  the  beam  element,  the  stiffness  matrix  is  four-by-four,  and  is  expressed  as 

"Kjfs)  K^fs)  -Kj($)  K3(5)- 
K<(s)  K2(s)  -K3(s)  Ki(s) 

-Kjfs)  -Kjfs)  Ksfs)  -K^fs) 

.  Kjfs)  K,(s)  -K4(s)  K2(s)J 


where 

Ki(s)=^{sh-st) 
Kjfs)  =  ^  (ch  St  -  sh  ct) 
K3(s)=^-(ch-ct) 
K4(s)  =  ^  sh  St 
K5(s)=^(sh  +  st) 

K  j(s)  =  j  (ch  St  +  sh  ct) 
A(s)  =  ^(1-chct) 


(2.66) 


(2.67a-g) 


This  sti'fness  matrix  also  reduces  to  the  static,  finite  element  stiffness  matrix  as  the  complex 
frequency  app.-oaches  zero.  The  ef.''ecis  of  disnibuted  forcing  and  initial  conditions  are  determined 
ria 


q<!(s)  = 


I 


r 

’a  g){x,s)‘ 

El 

■  E2(*'S) 

:a(I-chci) 

a  E3(x.s) 

V 

L  g4(x.s)J 

0 


f(j(x,s)dx 


(2.6S) 


and  the  Lnte.’polation  matrix  is  most  easily  expressed  by 


r(x.s)  =  0(x,s)['i'(s)]’'  - 


1x00 
0  10  0 
0  0  0  0 
LO  0  0  0. 


(2.69) 
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where 


and 


0(x,s) 


£«  e-“‘  e>“ 

ae®*  .ae-“‘  jae*"* 

EIo2e“*  EIa2e-“*  -Ela^ej"* 
.■Ela-’e“*  EIa3e-“*  jElaSt^* 


e'j®* 

-jae-j“ 

-Elo^e'j®* 

.jEJo3e-i»*. 


['i'(s)]’’ 


! 

api 

oL 

ae  % 

P2 

-e  p3 

opa 

/•.3L_ 
<xe  Pi 

P3 

,aL 
-e  P2 

4o(l  -  eh  Cl) 

aps 

-jP6 

oPS 

-jP7 

.oe^“‘'P8 

ae^“S5 

Lt  this  last  equation,  the  following  definitions  have  been  used: 

Pi  =  1  -  e'°^(ct-st)  P5  =  1  -  e‘j“^(ch-jsh) 

P2  =  !  -  e'°^(ci+si)  P6  =  1  -  e'j“''(ch+jsh) 

P3  =  e‘“^-(ct-st)  p7  =  e'j°^  -  (ch-jsh) 

P4  =  e'®^  -  (ct+st)  P8  =  e’j®'-  ■  (ch-jsh) 


(2.70) 


(2.71) 


(2.72) 


As  was  the  case  with  the  axial  rod,  die  expressions  for  the  internal  energy  of  deformaaon 
are  complex,  but  nonetheless  expressible  analytically.  The  potential  energy  maaix  is 


Hu(si.S2)  =  ai02^(si)‘’^ 


while  die  kinetic  energy  maaix  is 


FOl)  F(32)  -FOli)' 

?(•?:)  F(.3,)  -FOiSJ  .F(-?j) 

fO?4)  -FOPa)  F0?i)  FO?:) 
L-F(.jP3)  •F(.j?4)  Ff-I?,)  F(-j?,). 


V(S2)'' 


Ht(Si.S2)  =  'l'(Sl)’' 


-  4'(si)- 


F(Pi)  F(3j) 
F(-P:)  F(-?,) 
F0?a)  FOPj) 


LF(-j?3)  F(.j?4)  F(.j?2)  Ff.jpi), 


I 

ii  2  ii  3 
2^  ?- 
0  0 
LOO 


r 

G(Oi)  0  0' 

rFioj)  Ff-oj)  F(j02)  FC-jOj)-] 

F(.ai) 

G(-Oi)  0  0 

0(02)  0(-a2)  G(jo2)  G(-ja2) 

F(iai) 

GOoi)  0  0 

0  0  0  0 

LF(.jai) 

G(-ja,)  0  oJ 

L  0  0  0  oJ 

'J'(S2) 


(2.73) 


(2.74) 


where 
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(2.75) 


F{p)  = 

Om  = 

2nd 

P,  =aj  +«2 
p2  =  at  -  02 
P3  =  a, +j02  (2.76) 

p4  =  Oi  -  joi 

As  was  ihe  case  for  the  axiid  rod,  high  and  low  frequency  asymptotic  approximations  are 
required  for  the  beam  element  as  well.  These  approximations  allow  efficient  numerical 
computation  of  the  required  quantities  without  sacrificing  numerical  accuracy. 

2  5  6  Ktgh.frequency  elements 

Ln  many  cases,  the  simple  axial  rod  and  BeroouUi-Euler  beam  models  are  insufficient  to 
capture  die  high  frequency  bt'^avior  of  die  strucrural  system.  This  is  particularly  noticable  when 
analvair.g  the  wave- like  propagation  of  energy  associated  wmh  impulsive  disturbance  sources.  For 
these  situations,  more  accurate  strucrural  models  are  required. 

The  simple  axial  rod  model  presented  in  Section  2.3.4  predicts  dispersion.frce  propagation 
of  elasr.c  w  ax  es  at  all  frequencies.  For  rods  of  circular  cross  section,  a  bener  model  is  available  if 
a  radial  degree  of  freedom  is  introduced,  as  shown  in  Fig.  2-6.  This  is  the  basis  for  the  Mindlin- 
Hennr,ar.n  axial  rod  theory,  which  can  be  expressed  mathematically  by  the  following  system: 

-  a7().+2G)^  vi(x,s)  +  pa7s7vi(s,s)  =  2aX^V2(x,s)  (2.77a) 

“  .laK^x|^  V)(x,s)  (2  77b) 

Li  diese  equations,  vj  and  V2  are  the  axial  and  radial  displacements,  respectively,  a  represents  the 
radius  of  the  rod,  X  and  G  are  the  Lame  constants  of  the  material,  and  k  and  Kj  are  empirically 
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Ux=V2(X,l)  Uj=-Vj(X,t} 


desennined  parameters.  This  model  )«Jds  the  dispersion  curves  shown  in  Fig.  2-7,  and 
reproduces  the  dispersion  chaiacterisrics  of  the  rod  more  accurately  than  the  simple  rod  model. 

The  modified  stiffness  matrix  for  the  Mindlin-Henmann  mv  del  is  given  in  Appendix  A. 

A  more  accurate  beam  model  is  the  Timoshenko  beam,  shown  in  Fig.  2-8.  This  model 
allows  for  shearing  of  the  cross  sections  with  respect  to  the  center-line  of  the  beam,  and  accounts 
for  rotaiy  inertia.  The  equation  of  motion  is  given  by 

S^v(x.t)  +  pA<^x,t)  -  pl[l+f]0‘<<x.O  +  ^v(x,t)  =  f<,(x.t)  (2.78) 

where  G  is  the  shear  modulus  of  the  material  and  k  is  an  empirically  determined  correction  factor. 
The  Timoshenko  model  is  capable  of  supporting  both  a  shear  and  a  bending  mode  of  propagation, 
as  shown  in  the  dispersion  curves  in  Fig.  2-9.  This  model  also  places  a  finite  upper  limit  on  the 
fiexural  propagation  speed,  which  is  mireaUstically  unbounded  in  the  BemouUi-Euler  model. 
Details  on  die  stiffness  matrix  for  the  Timoshenko  model  are  presented  in  Appendix  A. 

2.4  Two-dimensional  elements 

Two  dimensional  elements  a.-e  raodeUed  using  partial  differentia]  equations  with  diree 
independant  variables  (two  spatial  dimensions  and  time).  Therefore,  the  modelling  of  two- 
dimension  elements  using  the  TE.M  methodology  cannot  produce  exact  results,  as  was  the  case  for 
one-dimensional  models.  The  reason  is  that  the  djuairn'cs  equation  remains  a  partial  differential 
equation  in  two  spatial  variables  after  the  Laplace  transform  operation.  As  a  result,  there  exist  an 
infinity  of  homogeneous  solutions  for  any  particular  element  model.  To  this  infinite  set,  these 
corresponds  an  infinity  of  points  along  the  boundary  (which  spans  a  one-dimensional  domain)  on 
which  boundary  conditions  must  be  satisfied.  Nevertheless,  by  using  a  sufficiently  large  number 
of  homogeneous  solutions  and  considering  only  a  sufficiently  large,  but  finite,  set  of  boundary 
points,  an  accurate  TE.M  solution  is  possible.  An  example  of  boundary  discretization,  for  the  case 
of  in-plane  deformation,  is  shown  in  Fig.  2-’0.  It  is  conjectured  that,  for  a  given  amount  of 
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shear 


Pig.  2-9:  Timoslicr.ko  beam  di.spcrsion  cliaraclcrislics  (a|=4xl0'^,  ai.=2.8):  (a)  Dispersion 
eurves,  (b)  Pretpiency  spcelnim. 


(x3.y3)  (x6,y6) 


Fig.  2-10:  Topical  two-dknensionaJ  elemem  wiib  three  generalized  displacements  and  forces  at 
each  boundary  point.  In-plane  deformation  is  assumed.  The  numbering  .scheme  shown  continues 
along  the  entire  boundary  of  the  element. 


35 


r 


‘  P  aona]  capabiliiy,  dijj  approach  jieJds  rciulK  saperfor  in  accuracy  lo  ihe  nadinonal  FEM 

roeihodoJo^y.  We  consider  wo  two-dimensional  elements  here:  a  plate  bending  element  and  a 
plane  stress  eleiaent 

2.4.1  Plate  bending  element 

Tht  TEM  fomtularion  for  plate  elements  bears  some  resemblance  to  presious  work  by 
Kulla  (19P0).  The  transfonned  erjuation  of  motion  for  a  plate  in  bending  is 

DV^’(x,y,s)  +  ms2v(x,y,s)  =  0  (2.79) 

where  v  is  the  denecaon  no.-mal  to  the  plane  of  the  element,  m  is  the  mass  per  unit  area,  and  is 
the  biharmonic  operator.  Also.  D  is  the  bending  rigidity,  given  by 


12(1.v2) 


(2.S0) 


'^here  h  is  the  thickness  of  the  plate  and  v  is  the  Poisson's  ratio  of  the  materia).  .Vote  that  we  hase 
assumed  no  dismbuted  fo..cing  and  aero  initial  conditions.  The  homogeneous  solution  vector  ,s  of 
Infinite  dimension,  with  each  ennn’  hating  the  form 


Substituting  this  erjuation  into  Eq.  (2.79)  tields  the  characteri.stic  equation 


[«*  P?] 


=  0 


(2  SI) 


(2.S2) 


Thus,  for  each  compiex-valued  a  (or  P).  dtcre  exist  four  independent  homogeneous  solutions 
corresponding  to  tite  four  compiex-valued  P's  (or  a's)  obtained  fmm  the  characteristic  equation. 

Toobtar'n  an  approximate  plate  solution,  we  must  selecafinite  set  of  values  foraCorP). 

A  general  method  of  selecting  this  set  has  not  been  developed  in  this  research.  From  this  set.  we 
obtain  a  truncated  solution  vector,  and  the  approximate  expression  for  the  deflection  field  becomes 

v(x.y,s)  =  v„(x.y.s)^a(s)  (,  g3) 
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The  problem  has  now  been  reduced  to  determining  the  coefficient  vector,  a{s).  in  terms  of  the 

boundary  conditions.  A  finite  set  of  boundary  points  on  the  element  is  therefore  selected.  The 

total  number  of  boundary  conditions  specified  at  these  points,  n,  must  equal  the  dimension  of  the 
homogeneous  solution  vector,  so  that  the  problem  is  not  underspecified.  three 

boundary  conditions  are  imposed  at  each  point.)  The  boundary  condition  constraints  can  be 

wTinen  m  the  form 

W.(S)  =  LS[v(xpy,s)]  =  I®[vH(Xi,yi.s)f  a(s).  i=l "  (2.S4a) 

o,(s)  =  D«[v(Xi,y,s)]  =  D®[v„(Xi,yi,s)]^(s)  ,  i=l a  (2.84b) 

where  is  a  linear  spatial  differential  operator  (independent  of  s)  relating  the  approxunate 
solution,  V,  to  the  i'th  generalized  displacement  on  the  boundary,  wj.  Similarly.  'D^  relates  v  to 
the  corresponding  dual  generalized  force,  o^.  Grouping  Eqs.  (2.84a)  and  (2.84b)  into  matrix  form 

yields 

w(s)  =  'l'(s)a(s),  q(s)  =  '?(s)a(s)  (2.8oa,b) 

where 

l4,y  [''h(x  1  -yi.s)]^  D'J;y[vH(x  1  .y  i.s)] 

^  >5(5)  4  ■  (2.S6a,b: 

Finally,  the  d>-namic  stiffness  matrix  follows  from  eliminating  a  from  Eqs.  (2.8Sa)  and  (2.85b)'. 

K(s)  =  y(s)['i'(s)]'’  (2,87 

Obviously,  the  choice  of  boundary  points  affects  the  accuracy  of  the  solution,  and  should 
depend  on  the  geometry  of  the  element.  Unfonunately,  a  quantitative  relationship  between  the 
boundary  point  locations  and  the  solution  accuracy  has  not  yet  been  developed.  However,  a 
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general  rule  of  thumb  is  to  space  the  boundaiy  points  more  closely  together  near  comers  of  the 
element,  where  internal  stress  gradients  are  large. 

2.4.2  Plane  stress  element 

^°^P’^**««P™Mem,theclifficuItyHesinexpressingthesta.eofd^^^^ 

of  a  single  scalar  function,  as  there  arc  two  in-plane  deformanonal  degrees  of  freedom.  To  remedy 
.his  simadon.  it  becomes  nece  t  .y  ,o  derive  a  frequency-dependant  stress  funcuon,  from  which  all 
internal  stresses  and  both  displacements  can  be  determined.  This  is  accomplished  by  working 

&cctlyuith  the  basic  plane  stressrelations.  expressed  in  he  frequency.domain.F^ 

ijowpic.  plane  stress  p.'oblems.  the  smess-strain  relations  are 


^  =  i(<^X-VCy) 

(2  88a) 

£y  =  g(<^y-VCr,) 

(2  SSb) 

P  ->  -L-r 

^y  "  20  ^»y  ~  “V 

(2.S8C) 

-here  E,G,andvare  he  extensional  modulus.  Shear  modulus  andPoisson'srati-ooftbema,^ 

«-‘pecnt.ely.  The  equarions  of  force  equilibrium,  expressed  in  he  frequen'cy-domain.  are 
•af  +  ^  =  ps2u,.  +  ^  = 

he.e  p  IS  the  density  of  the  materiai  and  s  is  the  (generally  complex)  Laplace  vznsbh.  To  these, 
wc  must  add  the  geometric  relations 


c 


2nd  the  compatibility  constraint 

3x1  -  2-^ 

Substituting  Eqs.  (2.g8a)  and  (2  SSb)  into  equation  (2.91)  yields 


^  ay’  -  2(35?+  (2.90a.b.c) 
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3y2  ^  3*2  '^1,3x2  gyt)  3x3y 

while  adding  the  derivative  of  Eq.  {2.89a)  with  respect  to  x  to  the  derivative  of  equation  (2.89b) 


with  respect  to  y  produces 


Eiiminating  tj;y  between  the  previous  two  relarions  results  in 

Substituting  Eqs.  (2.88a).  (2.88b).  (2.90a)  and  (2.90b)  into  equation  (2.94)  yields 

+  +  ^  +  (l+v)(l-v)^(ax+  Jy) 

3)2  3x2  3x2  3yi  fi 


which  reduces  to 


[v2-(l-v2)e^]  (Ox  +  Oy)  =  0 


where  V2  is  the  two  dimensional  spatial  Laplacian  operator.  Because  C,  and  Oy  are  independent 
variables,  we  sttll  need  one  additional  equation,  so  as  to  uniquely  identify  Ox  2«<3  Th'  second 
equation  is  obtained  by  subtracting  the  derivative  of  Eq.  (2.89b)  with  respect  to  y  from  the 
derivative  of  Eq.  (2.89a)  with  respect  to  x,  yielding 


ps2  /  V 

(Uv)%(crx-<ry) 


This  relation  is  equivalent  to 


2/  X  •)/  \  82ox  g2q^  Shy 

2(1+V)f  -  (Ox  '  Oy)  =  C<2x  3y2  +  3x2  '  ay2 


from  which  we  obtain 


[V2  -  2(1 +v)  ef]  (Ox  •  <2y)  =  [-  §  +  ^](0x  +  ^y) 
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Equations  (2.96)  and  (199)  represent  two  linear  partial  differential  equations  in  two 
variables  (the  sum  and  difference  of  the  normal  stresses).  This  system  can  be  reduced  to  a  single 
equation  by  defining  a  frequency  dependent  stress  function,  <!>,  so  that  the  following  relations  hold: 


(CX  +  Oy)  =  [v2  -  2(l+v)  O 

-^y)  =[- 


(2.100a) 


(2.100b) 


Then,  in  terms  of  <I>,  Eq.  (2.99)  is  identically  satisfied,  while  Eq.  (2.96)  becomes 


iIa'  -  Vf-  '^0  *  VF 


(2.102a,b,c) 


The  constants  Cj,  Cj  and  Cq  are  readily  identified  as  the  propagation  velocitie.  of  compression  and 
shear  wa\’es  in  a  plane  and  compression  waves  in  a  three  dimensional  medium,  respectively.  It  is 
interesting  to  note  that,  under  static  conditions  (s=0),  equation  (2.101 )  reduces  to  the  familiar 
biharmonic  equation  associated  with  the  Airy  stress  function. 

It  remains  to  determine  the  physical  entities  of  interest  in  terms  of  O.  In  the  traditional 
stress  function  formulation,  the  stresses  are  expressed  as  derivative  operators  on  O.  In  order  to 
obtain  similar  differential  operator  expressions  in  this  development,  it  is  necessary  to  define  a 
related  function,  v,  such  that  the  following  relationship  holds: 


The  final  form  of  Eq.  (2,101)  then  becomes 


Equations  (2.100)  and  (2.103)  can  then  be  used  to  determine  the  normal  stresses  in  terms  of  v. 
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=  =  (2.1053, b) 

Now,  combining  Eqs.  (2.88a),  (2.90a)  and  the  denvative  of  Eq.  (2.89a)  with  respect  to  x  leads  to 

=  f  (0,-vay)-^  (2.106) 


Subsntuung  Eqs.  (2.105a)  and  (2.105b)  into  equation  (2.106)  and  integrating  with  respect  to  both 
X  and  V  yields 


=  r  _3L.  +  dv2.^1, 

[  ox^ayi  c§  c3c|J 


Thus,  Eqs.  (2.105)  and  (2.107)  are  the  des'red  relations  between  the  stresses  and  v.  Furthermore, 
making  use  of  Eqs.  (2.88)  and  (2  90)  leads  to 


ax2  ciJ'' 

(2.108a) 

ii.dlv 

clJ'' 

(2.108b) 

Thus,  unUe  the  traditional  stress  function  formulation,  this  development  also  expresses  the 
displacements  in  terms  of  the  frequency  dependent  funcdon,  v. 

^\'e  are  now  in  a  posidon  to  apply  the  same  methodology  as  was  used  for  the  plate  bending 
clement.  Once  again,  the  assumed  form  of  the  homogeneous  solution  is 

VH,(x.y,s)  =  i=l,...,n  (2.109) 

where  a  and  p  are  funcaons  of  s,  leads  to  the  characteristic  equation; 

«iPi[(a|2+P,^)  -  ^j|^(a,2+p,2)  -  =  0,  i=l,...,n  (2.110) 


Of  the  four  solutions  to  this  equation,  two  (u=0  and  p=0)  are  spurious.  The  other  tw  o  determine 
the  relationship  that  must  hold  between  a  and  p  for  each  basis  solution.  The  rigid  body  modes  are 
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accounied  for  by  setting  s  equal  to  zero.  As  is  the  case  for  the  plate  bending  element,  the  acmal 
choices  for  ct,  and  P,  vary  depending  on  the  geometry  of  the  element  and  are  not  discussed  here. 


2.5  Assembly  Procedure 

The  assembly  of  indi\’idual  element  to  fonn  a  global  structural  model  is  performed  in  a 
manner  similar  to  traditional  FEM  techniques.  The  element  equations  are  first  collected  into  a 
large,  unreduced  matrix  equation,  given  by 


q‘(s) 

'k’(s)  , 

0  ■ 

■«'($)■ 

'qj(s)' 

= 

K^(s) 

+ 

-q^’(s)- 

.  0 

K'‘'(S). 

-w>’(s)- 

1 - 

JO 

w  here  the  superscripts  identify  the  indi\’idual  eleme.nts.  The  geometry  of  die  interconnecrions 
between  elements  is  specified  by  a  connectivity  matrix,  C,  which  relates  the  local  boundary 
displacements  of  the  elements  to  a  set  of  global  displacements,  w°(s),  diat  define  the  global  model. 


'w’(s)' 

.«>'(s)- 


=  CwG(s) 


(2.112) 


Figure  2-11  presents  a  simple  structural  system  and  the  associated  connectivity  matrix.  Piche 
(19S63)  shows  that,  for  small,  linear  displacements,  die  following  dual  relationship  holds: 


q?(s)  = 


q’(s)' 


Lq-V 


(2.113) 


Using  these  connectivity  relations  in  Eq.  (2.1 1 1)  yields  the  unreduced  system  model 

q°(s)  =  K^(s)vvC(s).^q°(s)  (2.114) 

where 
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Fig.  2-11:  Assembly  of  simple,  four  element  stmerure.  Elements  1  and  3  are  axt'aj  moocls.  and 
ele.ments  2  and  4  are  beam  models.  The  joint  is  assumed  to  be  massless. 
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(2. II 5a, b) 


K^(s)  = 


‘k’(s)  , 

0  ■ 

k\s) 

c. 

q  =  c"" 

.  0 

At  this  point,  if  there  are  any  global  displacements  that  are  constrained,  they  are  removed  from  the 
unreduced  model.  This  produees  the  desired  global  model,  given  by 

q%)  =  K°(s)  w°(s)  +  q°(s)  (2.116) 

Equation  (2.1 16)  represents  the  global  dynamic  stiffness  matrix  for  the  structural  model.  It 
is  mathematically  exact,  and  must  be  calculated  at  each  frequency  of  interest.  To  compute  the 
response  of  the  model  to  various  excitations,  we  must  solve  for  the  global  displacements.  Thus, 

=  G(s)[q°(s)-q5(s)].  G(s)  =  [K°(s)]’'  (2.117a,b) 

Here,  G(s)  represents  the  global  transfer  function  matrix  for  the  model.  If,  in  addtrion,  the  local 
boundaiy  displacements  for  a  particular  ele.ment  are  desired,  a  partition  of  the  connecti\'ity  matrix 
must  be  used: 

wks)  =  C^G(s)[q°(s)-q°(s)]  (2.118) 

FLnally,  the  internal  states  of  a  particular  element  a.'e  available  xia 

u>(x,s)  =  r(x,s)C4G{s)[q°(s)-q°(s)]  +  u’(x.s)  (2.119) 

bi  practice,  the  matrix  multiplications  are  never  performed  literally.  The  conneciixity  and 
unreduced  stiffness  matrices  are  highly  structured,  malting  it  possible  to  write  specialiaed 
algorithms  for  each  equation  given  abox'c.  This  dramatically  increases  the  overall  computation 
speed  of  the  assembly  process. 


2.6  Applications 

This  secrio.n  discusses  some  of  the  applications  of  the  TEM  modelling  approach.  Ail 
applicanons  ictjuire  global  model  assembly,  and,  consequently,  a  genetal  TEM  structural 
modelling  code  has  been  developed. 

2.6.1  Modal  frequencies 

In  many  cases,  all  that  is  required  from  a  structural  model  is  a  set  of  modal  frequencies. 
While  the  TEM  methodology  provides  significantly  more  information,  it  is  nonetheless  poss.ble  to 
obtain  modal  frequencies  using  an  algorithm  developed  by  Wittrick  (1971).  This  robust  algorithm 
uses  info.'mation  about  the  stiffness  matrix,  evaluated  at  a  trial  frequency,  to  determine  the  nu.mber 
of  modes  whose  frequencies  are  below  the  trial  frequency.  .Also  required  by  the  algorithm  are  uhe 
modal  frequencies  of  the  indixidual  elements  with  all  boundary  displacements  consaained  to  zero. 
The  algorithm  is  designed  for  undamped  structures  only,  and  additional  root  searching  techniques 
are  requL'ed  Ln  the  analysis  of  damped  structures.  Even  in  damped  cases,  howex’er,  the  algorithm 
protides  a  reasonable  initial  esdmate  of  the  location  of  the  damped  modal  frequencies.  Tnis 
algorithm  was  not  implemented  in  this  research,  as  adequately  accurate  modal  information  was 
available  from  plots  of  appropriate  transfer  functions,  as  described  below. 

2.6.2  Frequency  response  and  transfer  functions 

Tne  primary  advantage  of  the  TE.M  approach  is  its  ability  to  provide  the  exact  transfer 
function  matrix  at  any  frequency  of  inte.'-esL  This  is  obtamed  by  numerically  inverting  the  dt’namic 
stiffness  matrix.  The  stiffness  matrix,  K(s),  represents  a  matrix  of  coi.'.plex  impedances  relating 
generalized  boundary  forces  to  boundary  displacements.  Consequently,  G(s)  is  a  matrix  of 
complex  adminances,  and  is  often  called  the  dynamic  flexibility  matrix. 

The  transfer  functions  of  cantileveled  axial  rods  and  Bemoulli-Euler  beams  with  various 
damping  models  are  shown  in  Fig.  2-12.  For  the  rod.  the  input  is  an  applied  force  rn  the  free  end 
and  ihe  output  is  the  axial  deflection  at  that  end.  Similarly,  for  the  beam,  the  input  is  an  applied 
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(a) 


damping  (-  •  -).  and  wth  square  root  damping  (-  -  -):  (a)  Axia)  tip  force  to  axial  tip  displacement, 
(b)  Transverse  tip  force  to  transverse  tip  displacement. 


transverse  force  on  the  free  end  and  the  output  is  the  transverse  deflection  at  that  end.  It  should  be 
noted  that  the  accuracy  of  these  transfer  functions  extends  arbitrarily  high  in  frequency,  insofar  as 
the  inathetnatical  models  represent  the  actual  physical  system. 

In  order  to  demonstrate  the  capabilities  and  advantages  of  the  TEM  methodology,  a 
reasonably  complex  strucnire  was  anal)^^.  The  Spacecraft  COntrol  Laboratory 
Experiment  (SCOLE)  model  is  a  three  dimensional  asymmmetric  structure  proposed  by 
NASA  as  a  design  challenge  by  Taylor  (1986).  It  consists  of  a  rigid  shuttle  and  hexagonal 
truss  antenna  connected  by  a  flexible  mast,  as  shown  in  Fig.  2-13.  PreWous  authors  have 
treated  the  antenna  as  being  rigid.  In  this  effort,  however,  the  flexibility  of  the  antenna  is 
considered.  The  TE.M  model  thus  contains  thineen  beam  elements  (the  mast  and  nvelve 
antenna  elements)  and  a  six  degree  of  freedom  rigid  mass  representing  the  shuttle.  In 
addition  to  the  six  rigid  degrees  of  freedom,  a  total  of  52  partial  differential  equations, 
incorporating  axial,  bending  and  torsional  modes,  are  modeled.  For  comparison,  the 
SCOLE  model  was  also  analysed  using  ASTROS,  which  incorporates  a  finite  element 
algorith,m  similar  to  that  found  in  N.ASTRAN.  For  the  finite  element  model,  the  mast  uas 
divided  into  32  equal  elements,  and  each  of  the  antenna  beams  was  divided  into  four 
lumped  elements,  leading  to  480  degrees  of  freedom. 

Figure  2-14  compares  the  mansfer  functions  from  a  torque  applied  to  the  shurJe  about  the 
axis  of  the  mast  to  various  points  along  die  mast  and  antenna.  The  TEM  and  FEM  models  agree 
rather  well  at  low  frequencies.  However,  it  is  clear  that  die  finite  element  model  becomes 
inaccurate  beyond  the  first  few  modes.  WTiat  is  considerably  more  striking  is  the  relative 
computation  time  requred  to  generate  the  transfer  functions  shown.  On  a  micro- VAX  machine,  die 
TEM  analysis  required  approximately  one  hour  of  CPU  time,  m  contrast  with  several  days  of  CPU 
time  for  the  FEM  approach.  This  remarkable  acceleration  is  due  primarily  to  the  reduction  of  the 
total  degrees  of  freedom  in  the  model,  which  is  associated  with  the  lack  of  spatial  discretization  of 
the  beam  elements.  Since  the  compulation  time  associated  with  matrix  inversion  is  roughly 
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Fig.  2-13:  Scliciiialic  of  llic  SCOLli  .slrucUirjl  system  will)  flexible  nnleiiiin  model.  Node  numbers 
are  indicated. 


proportional  lo  the  cube  of  the  dimension  of  the  matrix,  the  reduction  in  total  degrees  of  freedom 
has  a  profound  effect  on  coroputaion  time. 

2.6.3  Time-domain  simulation 

A  final  application  of  the  TEM  approach  for  structural  analysis  is  in  the  time-domain 
simulation  of  structural  responses.  This  is  accomplished  via  the  inverse  Laplace  transfonn 
algorithm  presented  in  Sec.  2.2.  The  flexibility  matrix  is  evaluated  at  a  finite  set  of  N  frequencies, 
and  is  multiplied  by  the  global  force  vector,  which  contauis  the  Laplace  transforms  of  the  forcing 
functions  etaluated  at  those  same  frequencies.  The  resulting  displacement  vectors  are  then 
collected,  and  the  algorithm  generates  the  time-domain  responses  evaluated  at  a  set  of  2N  points  in 
time. 

rig.  2-15  compares  the  time-domain  simulations  of  a  simple  axial  rod  and  a  Mindlin- 
Hemnann  rod.  The  dispersite  effects  of  the  higher  order  model  are  apparent.  Likewise,  Fig.  2-16 
compares  the  Bemoulli-Euler  and  Timoshenko  beam  models.  Here,  the  effect  of  finite  disturbance 
propagation  \'elocity  is  tiie  primary  distinction  between  the  models. 
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3  MULTIBODY  TEM  FORMULATION 


The  TEM  foimularion  described  in  Chapier  2  pro\ndes  an  exact  PDE  mcxJel  of  complex 
sirucmres  with  small  motions.  Because  tiie  nodal  displacements  are  all  referenced  to  the  inertial 
frame,  the  small  deformation  assumption  for  the  elastic  deformation  also  implies  small  rigid  body 
morion  of  the  total  structure.  To  allow  larger  ranges  of  motion  at  articulated  joints,  one  can  embed 
the  element  reference  frames  at  the  undeformed  element  location,  which  may  have  both  rigid  body 
rotational  and  translational  motion.  This  U’pe  of  approach  has  been  used  in  multibody  tools  such 
as  DISCOS,  TREETOPS,  DADS  and  ADAMS.  As  shown  in  the  following  sections,  the  coupling 
of  the  rigid  and  elastic  degrees  of  freedom  results  in  a  set  of  integro-partial  differential  equations. 
The  fOi-mulation  will  be  derived  using  a  planar  example.  The  extension  to  three  dimensions  should 
be  straightforward. 


3.1  Mathematical  Model 

3.1.1  Equations  of  .Morion 

Consider  a  single  uniform  beam,  cantilevered  to  a  rigid  mass,  as  shown  in  Fig.  3-1.  The 
beam's  coo.-dinate  frame  is  fixed  to  the  rigid  mass,  which  can  undergo  rigid  body  motion.  For 
simplicity,  the  motion  of  the  system  is  assumed  to  be  planar.  Torques  and  forces  can  he  applied  to 
cither  end  of  the  beam'mass  system.  The  equations  of  motion  will  be  derived  using  Hamilton's 
principle. 

The  vector  from  the  inertial  origin  to  an  arbitrary  point  on  the  beam  is  given  by 

«■» 


where  all  uhe  vectors  are  expressed  in  body  coordinates.  The  velocity  of  the  point  is  given  by 


P  = 


L'yj 


[vy  +  a(x+v^)J 


(3.2) 


where 
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(3.3) 


The  kinetic  energy  for  both  the  rigid  mass  and  the  beam  can  then  be  written 


L 

T  =  |m[v2+v2]  +|mj  [V2  +  v2  +  2avyV,-  2v^avy 

+ Vy + Vy  +  + 2Vy Vy + 2o)(x+v^)Vy + 2Vya(x4v^)]  dx  (3.4) 


where  M  and  I  are  the  mass  and  inertia,  respectively,  of  the  lumped  mass,  and  m  is  the  mass  per 
unit  length  of  the  beam.  The  potential  energy  due  to  bending  and  axial  extension  is  given  by 


VL-tuai  work  due  to  the  external  loads  is  given  by 

5W  =  (F„-F,j)5X  +  (Fy,.Fj.j)5Y  +  [t,  ■  Tj +XFy,  •  V(F„.F,j)  -  (X+L)Fyj5e 

4  [yF,,-(X4L)Fyj-Tj]5v'(L}-Fy^6Vy(L)  (3.6) 

One  can  then  invoke  Harmlion's  priciple: 

jW-o>»w]d, .  0.  (3.7) 

Substituting  the  kinetic  and  potential  energies  into  Eq.  (3.7),  performing  the  appropriate  variational 
differentiation,  and  collecting  the  coefficients  of  5X,  5Y,  69, and  6vy,  one  can  obtain  the 
following  set  of  integro-partial  differential  equations: 

L 

(l4|mL5)a4|mL2Vy4m JxVy(x)dx  =  Ti-T2-LFy2  (3.8) 
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L 

(M+mL)Vx  +  mJvx(x)dx  =  FxfFxa 

(3.9) 

L 

A 

|^)L2o)•^(M+mL)Vy-^-m J  Vy(x)dx  =  Hyj-Fyj 

(3.10) 

m  Vx  +  m  v^(x)  =  EA  Vx(x) 

(3.11) 

*  *  •*  3^ 

mxco  +  mVy +  mvy(x)  =  -EI^^Vy(x) 

(3.12) 

The  nonlinear  terms  have  been  ignored  in  the  above  equations,  assuming  small  rigid  body 
velocities  and  small  elastic  deformations. 

3.1.2  Solution  for  the  Integral-Partial  Differential  Equations 

Equations  (3.S)  through  (3.10)  are  integral-differential  equations,  and  Eqs.  (3.11)  and 
(3.12)  are  partial  differential  equations.  The  set  of  equations  can  be  reduced  by  transforming  to  the 
Laplace  frequency-domain,  solving  Eqs.  (3.11)  and  (3.12)  forv^  and  Vy  in  terms  of  u,  V^,  and 
Vy,  and  then  substituting  the  result  into  Eqs.  (3.8)  through  (3.10). 

For  axial  extension,  let  us  define  the  axial  state  vector  in  Laplace  transform  space  by 

Ua(’‘)  =  [''^(x)  Ea|^\\(x)]^  (3.13) 

Equation  (3.1 1)  can  then  be  replaced  by  the  axial  state  equation 

|;Ua(x)  =  CaUJ(x)■^msVxP2  (3.14) 

which  has  the  following  matrix  exponential  solution: 

X 

U3(x,s)  =  e^2’‘ua(0,s)  +  msVx  P2  (3.15) 

where 
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Q(s)  =  [jjj  Pi  =  K]  .  P2  =  [?]  (3.16a.b,c) 

The  iitegial  of  v,.,  which  is  needed  in  the  transformed  version  of  Eq.  (3.9),  is  then  found  to  be 

LX  Lx 

Jvx(x,s)dx  =  pjje*'*’‘dx  UjCO.s)  +  p|msYjJJe^3^d4dx  P2  (3.17) 

Simiiariy,  for  transverse  bending,  let  us  define  the  bending  state  vector  in  transform  space  by 

Ub(x)  =  [vy(x)  |jvy(x)  El§vy(x)  El|^Vy(x)]^  (3.18) 


The  bending  equation  is  then  given  by 

|jUb(x)  =  Cb  Ub(x)  -  [m  X  s  <0  +  ra  sVy)  p4 

with  the  exponenriaJ  matrix  solution 


(3.19) 


X  X 

Ub(x,s)  =  Ub(0,s)  •  m  s  03  j  ^eCb(x-£)  sVy  J  d?  P4  (3.20) 


where 


Cb  = 


— 1 

0 

0 

0^ 

rii 

roi 

00^0 
0  0  0  1 

.  Pi  = 

0 

0 

.0. 

11 

0 

0 

.•ros^  0  0  0. 

(3.21a, b,c) 


The  two  integrals  of  Vy  in  the  Laplace  transform  of  Eqs.  (3.8)  and  (3.10)  can  then  be  wrinen  as 
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(3.22) 


LX  Lx 

I Vy(x,s)  dx  =  p]j dx  Ub(0,s)  -  p j  m  s  CJ  dx  P4 


Lx 

-  PjfflsVyJJe^I’^’^'^^d^dx  P4 


and 


L  X 


J xvy(x,s) dx  =  p[ J xe^l’’^ dx  U(,(0,s)  -  pjmscojxj d^ dx  P4 

L  X 

-  p[msVyJxJe^b(^'^)d5dx  P4  (3.23) 


The  Lniegral  expressions  are  substiiuted  into  Eqs.  (3.8)  through  (3.10),  and  the  coefficients 
of  SO),  sV^,  sVy,  s-Vj^,  and  sA’y  are  collected  into  a  matrix.  Equations  (3.15)  and  (3.20)  are 
evaluated  at  x=L,  the  left  and  right  sides  interchanged,  and  the  coe.'ficients  are  again  collected  into  a 
matrix,  O.  This  yields  an  equation  of  morion  of  the  form 


O 


s  y(s) 

s*Ua(0,s) 

.s2ub(0,s). 


=  f(s) 


(3.24) 


where 

y(s)  =  [o)(s)  V;,(s)  Vy(s)]’^  (3.25) 

f(s)  =  [Ti-TrLFyz  Fxi-Fxz  Fyj-Fy^  u.d^.s)”^  Ub(L,s)^^  (3.26) 

Let  us  now  replace  the  axial  states,  Uj(0,s),  and  the  bending  states,  Ub(0,s),  on  the  left  of 
Eq.  (3.24)  by  the  nodal  displacements  at  the  ro'o  ends  of  the  beam/mass  model.  The 


5S 


transforaatfon  is  done  by  pairirioned  matrix  solution  of  Eqs.  (3.15)  and  (3.20)  for  EAv^(0)  and 
[EIVy(0)  EIVy(0)],  respectively.  Asa  result,  the  axial  and  bending  partitions  of  the  foice  vector 
are  replaced  by  nodal  forces.  These  nodal  forces  are  equated  to  the  external  forces  and  torques  by 
the  equations 


(3.27) 


(3.28) 


Terms  in  Eqs.  (3.27)  and  (3.28)  involving  so,  sV^  and  sVy  are  moved  to  the  right  hand  side. 
Finally,  the  equation  of  morion  becomes 


O  sw  =  f 


(3.29) 


where 


tv=[/  s[v^(0)v^(L)]  s[vj.(0)  v;<0)  Vy(L)  v'(L))f 


fTal 

f  =H, 

■^Ho 

.‘yjJ 

1  1 

.P». 

H,  = 


Ho  = 


’1  0  0  0  0  -s2  0  0  01 

0  1  0-s2  0  0  0  0  0 

.0  0  1  0  0  0  s2  0  0. 

■-1  0  0  0  0  0  0  -s2  01 

0-1  0  0 -s2  0  0  0  0 
-L  0  -1  0  0  0  0  0  $2 


(3.30) 

(3.31) 

(3.32) 

(3.33) 
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3.1.3  Ojiuiecrion  :o  Adjacent  Elements 

Let  us  now  consider  the  beatn/mass  model  as  one  element  of  a  larger  truss/frame  structure, 
as  shown  in  Fig.  3-Z  The  external  forces  and  torques  at  the  two  ends  of  the  beam  can  be 
expressed  as  the  sum  of  applied  forces/torques  and  constraint  forces/iorques: 


Sij  + 


■Ta' 

Pxa 

R,, 


gOi  + 


(3.34) 


(3.35) 


where 


B 


Ij  = 


'10  O' 
0  cos  Bij  -sin  ejj 
.0  sin  0ij  cos  eij. 


(3.56) 


The  vector  g  contains  applied  forces/torques,  the  subscripts  (jj)  and  (qj)  indicate  inboard  and 
outboard  joints,  respectively,  of  the  j'lh  body,  the  X  vectors  are  the  constraint  forces/torques  at  the 
joints,  Bj,  is  a  matrix  of  kinematic  coefficients  which  transfoiros  the  constraint  forces/torques  from 
the  coordinates  of  the  Inboard  body  to  the  coordinates  of  the  current  body,  0jj  is  the  angle  at  the 
joint,  and  O  is  a  selection  matrix  which  picks  out  which  degrees  of  freedom  are  to  be  constrained  at 
the  joint.  Without  loss  of  generality,  Eq.  (3.36)  has  arbitrarily  assumed  that  the  constraint  forces  at 
the  inboard  joint  are  in  the  inboard  body's  coordinate  system.  Substitution  of  Eqs.  (3.34)  and 
(3.35)  into  Eq.  (3.29)  give,  for  the  j-th  body, 

OjSWj  =  gj  +  HijB^5),jXij  +  ]^Hoj<i>oj>^j  (3.37) 


where 


Sj  =  Hjjgjj  +  Hpjgoj 


(3.38) 
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The  summation  in  Eq.  (3.37)  is  cairied  out  over  all  of  the  outboard  Joints  for  tree  topologies.  Note 
that,  for  an  end  body,  there  is  only  an  inboard  hinge  and  no  outboard  hinge,  so  that  the  last  term  in 
Eq.  (3.37)  disappears. 

The  hinge  Idnematics  relate  the  absolute  velodties  of  adjacent  bodies  to  the  relative 
velocities  at  the  adjoining  hinge.  For  the  inboard  hinge  of  the  j-ih  body,  we  have 


[Boj  O]wj  +  Boj.,  [I  Pjwj.i 

(3.39) 

r  -1  0  O-i 

=  0-1  0 

(3.40) 

L-L  0-lJi 

rO  6  0  0  0  It 

010000 

Lo  0  0  0  1  oJ 

(3.41) 

The  matrix  Ci  is  a  selection  matrix  for  die  unconstrained  degrees  of  freedom,  ct  is  a  vector  of 
rheonomie  constraints,  P^is  a  vector  of  unconstrained  relative  hinge  velocities,  and  Boj  is  a  matrix 
of  Itinematie  coefficients  that  relate  body  velocities  to  the  outboard  joint.  For  the  base  body  (j=l), 
the  last  term  in  Eq.  (3.39)  disappean,  because  the  base  body's  inboard  body  is  the  inertial  origin, 
which,  by  definition,  has  aero  velocity. 


3.1.4  Recursion  Solution  for  the  Total  Structure 

We  now  combine  the  dynamics,  represented  by  Eq.  (3.37),  and  the  Itineroatics,  represented 
by  Eq.  (3.39).  Let  us  consider  an  end  body  (j=e)  and  solve  the  kinematics  equation  for  j. 

Ye  =  B,.-’  [<&i.  p,l  +  <i)i,6tfc  -  Bo.c.,  [I  P]  w,.,]  (3.42) 

Equation  (3.42)  is  augmented  by  die  nodal  velocities  and  re-wiinen  as 

rll 

We  =  'J'ePe  +  [0jB,;'  K  ctj.  -  Bo...,  [I  P]  w..,]  (3.43) 
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where 


'i'e  = 


B 


■h 


U  '^U  ®  ® 
0  I  0 

L  0  0  1. 


Pe  =  [Plf  s[Vjj(0)  Vj^(L)]  s[vy(0)  vj,(0)  Vy(L)  vJ.OL)]]’’ 


(3.44a) 

(3.44b) 


Equation  (3.43)  expresses  the  body  j  velocities  in  terms  of  the  velocities  of  its  inboard  joint  and 
inboard  body. 

Substinttion  of  Eq.  (3.43)  into  the  dynamics  equation,  pre-roultiplying  by  and  solnng 
for  Pj  gives  an  expression  in  terms  of  inboard  body  velocities  and  inboard  hinge  constraint 
forees/torques: 


rh 


[^!=aic-Bo.c.i[I  P]we.l]} 


(3.45) 


The  unknown  constraint  torques  can  be  solved  for  by  pre-multiplpng  Eq.  (3.39)  by  ij/  and 
substituting  Eq.  (3.37): 

h  =  [^i. Oj.  •  J  [Bic  0]  ie -  [I  P]  w,.,]  (3.46) 

where 

Je  =  [IK  [Bj.  0]  5;'  %  (3.47) 


The  unconstrained  relative  joint  velocity  term  has  disappeared  because  and  are  orthogonal 

to  each  other. 

Equations  (3.43),  (3.45)  and  (3.46)  are  the  recursion  relations  for  the  end  body.  Given  the 
velocities  of  body  c-1,  one  can  compute  the  constraint  forces/iorques  via  Eq.  (3.46),  substitute  into 
Eq.  (3.45)  to  get  the  relative  velocities,  and  substitute  into  Eq.  (3.43)  to  gel  the  body  velocities. 
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To  obtain  the  recursion  for  intermediate  bodies  in  the  middle  of  the  chain  or  tree  lopoiogy, 
let  us  consider  the  e-1  body  as  the  current  body  and  relate  its  velocities  to  its  inboard  Wnge.  The 
interconnection  constraint  forces/toniues  have  been  defined  so  that  =  X^j.  We  can  therefore 
substitute  Eq.  (3.46)  into Eq.  (3.37)  with  j  equal  toe-1,  and  rearrange  to  yield 

‘^e-l  S^e-l  =  St-l  +  ^i.e-i^^I.e-l  ^I.e-t  (3-^8) 

where 

<1  =  5e.l  +  X  s  [I  P]  (3.49) 

0 

gt'l  =  ie.l  +  X”0'->  -  7  fB,,  0)  O’J  ie)  (3-50) 

o 

Equation  (3.50)  is  now  in  the  same  form  as  for  an  end  body,  and  the  same  derivation  can  be 
followed  to  produce  a  recursion  for  the  e-1,  as  well  as  all  intermediate  bodies.  Recall  titat  at  the 
base  body  0'=1)  the  vector  wj.],  which  conresponds  to  the  inertial  frame,  disappears  because  tiie 
Ltertial  origiti  has  zero  velocity.  This  gives  the  initial  value  for  the  forward  recursion  to  begin  from 
the  base  body  to  the  end  bodies.  The  multibody  algorithm  thus  involves  aijacbv  ard  recursion  to 
compute  the  equivalent  matrices  and  vecio.-s,  O’  and  g’,  and  a  forward  recursion  to  compute  the 
Wj  vectors. 

3.2  Discussion 

The  mathematical  formulation  has  been  ptesenied  for  a  linearized  model.  The  linearization 
has  been  necessary  for  the  Laplace  transformation  to  be  applied  in  solving  the.  PDE  models.  It  is 
felt  that,  even  though  the  equations  have  been  linearized,  the  range  of  angular  motion  that  can  be 
simulated  has  been  enlarged  when  compared  to  the  cartesian-based  TEM  models.  The  derivation 
was  shown  for  a  planar  model  for  ease  of  presentation.  The  extension  to  iluee  dimensional  motion 
should  be  straightforward. 
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•Die  current  fonnuladon  is  applicaWe  to  chain  and  tree  topologies.  In  order  for  this 

forntulariontobeapph-ed  to  trusses  and  fiatnes.it  must  be  extended  tohandJecbsed 

loops.  This  tjpe  of  fomulation  has  been  performed  for  modal))-  based  dynamic  models  by  Chun 
(1991)  and  can  be  easily  adapted  to  PDE  models. 

Arbiranly  large  angularmodon  of  the  total  smicture.  as  well  as  amculated  joints,  requires  a 
nonlinear  model  for  dte  con-ect  description  of  both  the  dynamics  of  d,e  motion  and  the  kinematics 
at  thejoints.  Future  effons  should  explore  fte  use  of  perturbation  techniques  titat  allow  the 
Laplace  cransfomi  to  be  used  while  stiU  including  the  effects  of  the  nonlinear  tenns. 


4  CONTROL  DESIGN  BASED  ON  TEM  MODELS 

The  exactness  of  the  TEM  fomiulation  makes  it  possible  to  achieve  lemaikaWe  peribimance 
in  open-loop  slewing  maneuvers  of  flexible  structures.  This  is  the  subject  of  Section  4.1  of  this 
chapter.  Unfominately,  because  the  TEM  methodology  does  not  immediately  >16^  a  state-space 
representation  of  the  structual  system,  traditional  state-space  control  methods  are  not  directly 
appUcable  in  the  closed-loop  case.  (In  actuah'ty,  no  finite  representation  can  exist,  as  the  smicmral 
model  is  of  infinite  order.)  Methods  for  achieving  closed-loop  control  solutions  without  state- 
space  models  are  discussed  in  Section  4.2. 

4.1  Open-Loop  Control 

Ln  this  section,  we  develop  an  open-loop  control  algorithm  that  takes  advantage  of  the 
qualitj’  of  the  strucnsral  model  available  via  the  TEM  methodology.  We  restrict  anention  to  finite¬ 
time,  linear  maneuvers  with  a  quadratic  cost  functional.  We  also  assume  that  the  struemre  is 
initially  at  rest.  The  desired  terminal  state  is  expressed  by 

y(tf)  =  Td  (4-1) 

where  tf  is  the  maneuver  time,  y(t)  is  a  vector  of  variables  of  interest,  and  vj  contains  the  desired 
terminal  values  of  these  variables.  The  elements  of  y  could  include,  for  example,  the  displacement 
and  rotation  of  a  rigid  mass  on  the  structure,  or  the  relative  transverse  deflection  of  a  point  on  a 
flexible  member.  The  available  control  forces  are  also  collected  into  a  single  vector,  qc(t),  of 
dLmension  Nj.  These  control  forces  are  then  a  subset  of  the  global  generalized  forces  defined  by 
the  system  model.  In  order  that  these  forces  be  continuous  in  time,  we  must  impose  the  additional 
constraint 

qc(0)  =  <lc(0  =  0  (4.2) 

Also,  the  quadratic  cost  functional  is  given  by 
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I  =  I  [y(tf)-  yd]^R>-y[yOf)-  y<l3  +  5  J 

where  R„,  Rqq,  and  Rqq  are  s>-nimetric.  positive  definite  weighting  matrices. 

We  must  now  exFess  the  output  vector  in  terms  of  the  control  vector  vra  the  system 
dynamics.  This  is  accomplished  using  the  convolution  integral 

I 

y(t)  =  J  CyO-"')  qc(^) 

where  G/t)  represents  the  impulse  response  matrix  relating  y(t)  to  qj.i).  This  convolution  integral 
is  calculated  efficiently  using  the  inverse  Uplace  transform  algorithm  of  Sec.  2.2: 

y(0  -L-'[Gj.(s)qc(s)] 


4.1.1  BandOirruted  control  approximadon 

By  substituting  Eq.  (4.4)  in  Eq.  (4.3).  we  obser\-e  that  the  cost  functional  depends  on  q^ft) 
only.  Setting  the  first  variation  in  cost  to  aero  therefore  yields 


'  k 

f 

jGj.(tri)qc(')tit-.vaj  Ryy 

jGj.(trt)6qc(t)tii 

+  J[qc(')^Rqq5qc(«)  +  qc(0’'Rqq5qc(>)]d'  =  0 


Tlie  problem  then  lies  in  solving  this  equation  for  qc(t).  Unfortunately,  this  is  not  a  simple  maner, 
as  both  the  control  vector  and  its  variation  appear  within  the  integrals.  However,  if  the  control 
inputs  are  band-limited  (as  is  often  the  case),  a  numerical  solution  is  easUy  obtained.  Each  control 
input  is  first  approximated  by 

qci(0  =  fqW^Ci.  i  =  l . Nc 
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where  fq(t)  is  a  vector  of  known  basis  fiinctions  of  time  (usually  sine  and  cosine  functions),  and  q 
is  a  vector  of  undetennined  coefficients  coiresponding  to  the  i'th  control  input  The  entire  control 
vector  can  then  be  conveniently  expressed  as 

Ted)  =  Fq(t)c  (4.8) 


where  the  following  definitions  have  been  used: 


Fq(t)  = 


(4.9a.b) 


Using  this  band-limited  approximation,  the  optimal  control  problem  is  reduced  to  determining  the 
coefficient  vector,  c.  The  variations  in  the  control  vector  are  then 

5qc(0  =  Fq(05c,  Sq^CO  =  Fq(t)5c  (4.10a,b) 

and  the  constraints  given  by  Eq.  (4.2)  reduce  to 

Fq(0)c  =  Fq(tf)c  =  0  (4,11) 

Furthermore,  the  vector  of  desbed  outputs  is  expressed  by 

ydf)  =  V(tf)c  (4.12) 

where  Y(i)  is  the  basis  function  response  matrix,  and  is  given  by 

I 

Y(t)  =  jGj.(t-T)Fq(x)dt  =  L-'[Gy(s)F,(s)]  (4.13) 


4.1.2 


Solution  without  minimization 
Grouping  Eqs.  (4.1 1)  and  (4. 12)  yields  the 

rY(tf)l 

Fq(0)  c 

LFqdf)J 


matrix  equation 


Yd 

0 

0 


(4.14) 
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*  T 


ff  the  numb^rofdwimi  outputs  and  the  nuffiba  of  untaou,,  cocte 

tn  iii  ctjuadon  is  square,  then  c  is  uniquely  deteteuined.  Typically,  however,  there  art  many 

dements  in  c, so  tr,a.Eq.(4.i4)  is  underdetenntneAConseqt-en^^^ 

the  tenninal  constraints.  We  therefore  have  some  freedom  in  choosing  which  panicularc  to  use 

Foragiven  problem, the  part, cuiarchotcenumndres  some  predefined  co^^ 

pmWdas  a  measure  of  nominal  penonttanc  The  next  two  subsections  describe  -wo  such  cost 
funcdonals. 


4.3.3  ^^imizadon  with  point  consnainis 

Wefnst  use  theeost  fur, cfional  given  by  Eq.(4.3)and  adjoin  the  const^^^ 

(4.!  1)  wa  nvo  Lagrange  muldphers.  X,  and  X^  Taldng  variations  in  c  yields 
&}  =  [Y(tf)c  -  .v<i]’^RjyV(tf)Sc 

'f 

+  c  { J"  [PqW^RqqFqCt)  +  Fq(tl^R^(jFq(t)]  ot }  6c 

■^XjF,(0)5c.X]-F,(t^)6c^8XjF,(0)c  +  5X}F,(n)c  =  0 

leading  to  the  foUoving  matrix  equation: 

■  W  Fq(0)"  F(,(t,)“]r 1  fYr.vrp 

IfJ,,)  0  0  J[,;j  (  «  J  «.1«) 


tf 

■fhis  is  a  symmetric  system,  and  can  be  solved  using  standard  linear  algebra  routines. 

A  uPdque  advantage  of  this  approach  is  that  i,  readily  accomodates  penaioes  in  higher 
dertvauves  of  both  control  effon  and  structural  deformation.  In  the  frequency-domain. 
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(Ufferennation  merely  requires  multiplication  c  f  the  data  by  the  Laplace  transform  variable.  The 
inverse  transformation  then  produces  the  derivative  of  the  original  signal.  Higher  order  derivatives 
are  obtained  by  multipljing  by  higher  powers  of  the  complex  frequency.  Incoqxrrating  higher 
derivative  penalties  in  the  traditional  optimal  control  formulation  is  considerably  more  diffiailt. 

It  should  be  noted  that  the  only  approximatio.n  in  the  entire  development  involves 
expressing  the  control  inputs  in  terms  of  the  basis  functions.  The  dynamics  of  the  entire  structure 
are  accounted  for,  since  the  impulse  responses  are  exact  (insofar  as  the  original  equations  represent 
physical  reaiiw).  Also,  tire  structural  deformations  are  assumed  to  be  small,  so  that  linearization 
dees  not  introduce  significant  errors.  As  a  result,  large  angle  slew  maneuvers  are  not  included  in 
this  class  of  problems.  It  is  possible,  however,  to  express  structural  deformations  with  respect  to  a 
nominal  condition  during  a  large  angle  slew,  and  tlien  linearize  about  that  reference,  as  discussed  in 
the  presious  Chapter. 

In  an  earlier  analytical  snidy  by  Skaar  (1984),  the  open-loop  control  of  a  rigid  mass  with  a 
flexible  appendage,  shown  in  Fig.  4-1,  was  studied.  In  his  work,  structural  deformation  penalties 
were  not  incorporated  into  the  cost  function;  rather,  the  terminal  conditions  wtre  adjoined  to  the 
cost  functional  as  consaaints.  Skaar  deris  ed  analytical  expressions  for  impulse  responses  of  the 
simple  mass/appendage  saucture  and  thus  obtained  closed  form  optimal  conaol  solutions  for  the 
saucture.  Though  successful  for  this  application,  his  approach  does  not  readily  generalize  for 
more  complex  sauctures.  In  conaast,  the  formulation  presented  here  readily  generalizes  for 
realistic  complex  sauctures.  Skaar 's  example,  however,  is  used  as  a  first  example  to  validate  the 
optimal  conaol  formulation. 

The  maneuver  Lnvolves  translating  the  mass  a  distance  of  10  meters  along  the  axis  of  the 
flexible  appendage,  bringing  it  to  rest  with  minimal  residual  energy  and  post-maneuver  drift  after 
no  seconds.  The  first  case  places  terminal  penalties  on  the  final  position  and  velocity  of  tlie  rigid 
mass  and  on  a  point  4/5  of  the  length  along  the  flexible  appendage.  A  small  penalty  is  also  placed 
on  control  ra.e,  and  17  basis  functions  are  used  to  approximate  the  control  input.  The  results, 
shown  in  Fig.  4-2,  indicate  that  the  terminal  conditions  are  matched,  and  residual  energy  is 
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Ridd  Mass 


EA  =  0.05Kt 
EI  =  0.05Nt-m2 


Fig.  4-1:  Simple  mass/flexible  appendage  stnictural  model. 


negligible.  In  the  second  case,  the  member  sriffness  is  reduced  by  a  factor  of  four,  so  that  the 
piim^  modal  frequency  of  the  structure  corresponds  approximately  wth  thefttquency  of  the  first 
basts  function  of  the  control  inpuL  The  results  of  this  case,  presented  in  Fig.  4-3,  indicate  that  the 
control  input  has  been  adjusted  so  that  excitation  of  the  primary  mode  of  the  structure  is 
suppressed.  Again,  the  teiminal  conditions  are  matched,  and  residual  internal  energy  is  negligible. 

The  second  example  is  the  SCOLE  structure  analysed  in  Sec.  2.6.2.  The  maneuver 
presented  here  consists  of  a  ten  second.  0.1  radian  rotanon  about  the  z-axis  of  the  shutde.  Hiis 
maneuver  is  a  purely  academic  exercise,  and  is  unrelated  to  the  maneuver  specified  in  the  original 
design  challenge.  For  the  first  case,  torque  controls  directed  along  the  z-axis  are  placed  at  either 
end  of  die  mast  Due  to  the  asimmetty  of  the  structure,  gyroscopic  coupling  is  expected. 
Consequently,  roll  and  pitch  torque  controls  are  also  located  on  the  shuttle.  The  cost  of  control 
effon  is  equally  weighted  among  the  control  inputs.  Equal  terminal  magnitude  and  me  penalnes 
are  applied  to  the  roll,  pitch  and  yaw  angles  of  the  shuttle,  as  well  as  the  torsional  defoimadon  of 

the  2t  iis  siidpoint  and  at  the  msst/antenna  junction. 

The  results  of  the  first  SCOLE  slew  are  shown  in  Hg.  4-4.  It  is  clear  that,  although  the 
shuttle  has  rotated  the  prescribed  amount,  there  is  a  small  amount  of  residual  torsional  energy  m  the 
structure.  This  energy  is  due  primarily  to  the  deformation  of  the  antenna  and  mast  at  *e  lemunal 
time.  Also,  the  set  of  controls  utilized  are  L-.capable  of  suppressing  out-of-plane  deflecnon  of  the 

antenna,  which  is  caused  by  the  asymmetry  of  the  structure. 

In  order  to  suppress  this  residual  energy,  additional  controls  are  placed  on  the  antenna.  Iit- 
plane  femes  are  a  ’ailable  at  the  mast/an lenna  junction  and  directly  aaoss  the  antenna.  In  addition, 
an  out-of-plane  thruster  is  placed  at  the  laner  location.  Furthermore,  additional  penalties  are  placed 
on  antenna  deformation.  The  improvement  in  the  slew  response  can  be  seen  in  Fig.  4-5.  For  ,his 
maneuver,  most  of  the  torque  is  generated  by  the  antenna  thrusters  across  from  the  roast.  In 
reality,  this  distribution  of  control  effon  would  be  unwise,  as  it  would  lead  to  excessive  stress  in 
the  mast/amenna  junction.  Also,  as  shosvn  in  the  figure,  this  trajectory  causes  a  large  amotmt  of 
torsion? J  deformation  of  the  mast.  By  adjusting  the  relative  weights  on  the  controls  and  structural 
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deformation  outputs,  it  is  possible  to  converge  upon  a  more  realistic  trajectory.  However,  this 
control  solution  provides  an  adequate  demonstradon  of  the  formulation  presented  here. 

4.1.4  Minimization  of  flexible  energy 

Another  method  of  obtaining  an  optimal  solution  consists  of  minimizing  the  residual 
flexural  energy  within  the  structural  elements  at  the  terminal  time.  This  is  achieved  by  expressing 
the  generalized  boundary  displacements  of  the  i’lh  element  in  terms  of  the  undermined  coefficients: 

wi(s)  =  Ci  G,(s)  Fq(s)  c  =  Hi(s)  c  (4.18) 

Making  use  of  Eq.  (2.43)  then  yields 

E>‘(t)  =  |c'^Ei(t)c  (4.19) 

where 

09 

Ei(t)  =  k(t)  J  J Hksi)’’  Sksi.sa)  Hks2)  eJ<“i*“2>‘da>idc02  (4.20) 


Included  m  the  cost  functional  are  the  weighted  penalties  on  residual  energy  for  a  set  of  Xf  flexible 
elements  and  weighted  penalties  on  control  effort  and  control  rate.  To  this  we  adjoin  the  desired 
tentunal  conditions  and  the  constraints  on  the  controls  at  the  begirjiLng  and  end  of  the  maneuver. 
The  cost  functional  is  thus 

J  =  SriEktf)  +  ^^[qcCO^RqqqcCO+qcW^RqqqcWjdt 

+  ^'^[y(tf)-.v<l]  +  ^Jqc(0)  +  (4.21) 

Setting  variations  in  J  due  to  c  to  zero  yields 
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(4.22) 


■  W  Y(t/Fq(0)TF,(tf)'^ 

“c" 

X 

"O' 

Y(tf)  0  0  0 

yd 

Fq(0)  0  0  0 

0 

LFqCtf)  0  0  0  J 

L^fJ 

.0. 

where 

Kf  /•' p  .. 

W  =  X  n  +  J  L  FqCO^RqqFqC)  +  Fq(t)’‘RqqFq(t)J  dt  (4.23) 

i=l  ° 

Again,  this  system  is  symmetric,  and  can  be  solved  with  standard  linear  algebra  software 
packages. 

Tlie  imnimum  residual  energy  approach  was  applied  to  the  simple  mass/appendage  system 
studied  in  the  previous  section.  Two  maneuvers  were  performed,  both  with  a  prescribed  final 
displacement  of  10  meters  after  20  seconds.  In  the  first  maneuver,  the  desired  final  velocity  of  the 
rigid  mass  was  zero,  while  in  the  second,  the  final  velocity  was  1  meter/second.  Tne  results  of 
these  maneuvers  are  shown  in  Fig.'s  4-6  and  4-7.  In  each  case,  the  residual  energy  is  seen  to  be 
negligible. 

The  same  strucrare  was  used  to  perform  rotational  maneuvers.  In  Ms  case,  the  bending  of 
the  flexible  appendage  was  considered.  A  0.1  radian  slew  with  both  zero  and  1  radian/second 
terminal  angular  velocity  were  studied.  The  results,  shown  in  Fig.'s  4-8  and  4-9,  indicate  that 
perfonnance  comparable  to  the  axial  cases  was  achieved. 

The  minimum  energy  cost  functional  leads  to  system  trajectories  with  far  less  residual 
energy  than  those  obtained  via  point  constraints.  Funhermore,  minimization  of  total  deformational 
energy  also  avoids  the  problem  of  selecting  which  points  to  constrain.  All  that  is  required  is  a 
relati\-e  cost  weighting  for  each  flexible  element  of  interest.  However,  because  the  calculation  of 
inssmal  energy  involves  a  double  integral,  the  minimum  energy  approach  requires  more 
computational  effort  The  minimum  energy  cost  functional  was  not  applied  to  the  SCOLE 
maneuver  problem. 
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r ig.  4-6:  Linear  slew  maneuver  with  residual  energy  cost  runctional,  Icnninal  time  of  20  see, 
icnninal  displacement  of  10  m  ami  terminal  velocity  of  0  m/scc.  Control  forec  was  approximated 

with  18  basis  functions:  (a)  Position  of  rigid  mass  ( - )  and  tip  of  flexible  appendage  ( — ),  (b) 

control  force  applied  to  rigid  mass,  (c)  axial  deformation  of  flexible  appendage  at  center-span  ( — 
and  tip  (-  •  -)  with  respect  to  position  of  rigid  mass,  (d)  same  plot  as  (c)  with  expanded  vertical 


Displacement  (m) 


Displacement  (m) 


4.2  Closed-loop  control 

The  closed-loop  control  of  infinite  csder  systems  expressed  in  the  transformed  domain  is  a 
considerably  more  difficult  ptpblem.  The  exact  dynamics  are  available  in  die  fiequency-domain 
only,  and  no  finite  dimensional  state  space  realization  is  possible.  Asa  result,  full  order 
techniques,  such  as  LQR  and  LQG  methodologies,  are  not  applicable.  However,  the  control 
problem  can  be  posed  in  a  form  amenable  to  fiequency-domain  design  techniques.  It  is  assumed 
that,  for  a  given  smicnnal  model,  a  set  of  disturbance  forces  act  at  global  element  junctions,  and 
performance  is  measured  in  terms  of  some  set  of  global  generalized  displacements.  The  control 
objecdve  is  then  to  minimize,  in  some  sense,  the  transfer  funcdon  fi-om  the  set  of  disturbances, 
w(t),  to  the  performance  measure,  z(t).  This  is  to  be  accomplished  by  a  finite  order  controller 
which  has  available  as  inputs  a  finite  set  of  measured  generalized  displacements,  y(t),  and  acts  on  a 
finite  set  of  acmators,  u(i),  located  on  the  strucmre.  The  situadon  is  depicted  in  Fig.  4-10.  The 
transfer  funcdons  from  disturbances  and  control  inputs  to  the  performance  meoric  and  measured 
outputs  are  easily  obtained  as  parddons  of  the  dynarm'c  flexibility  matrix.  Note  that  this  control 
problem  is  in  the  "standard  form,"  which  has  been  studied  extensively  by  Francis  (1987)  and 
Doyle  (1989)  for  finite  dimensional  plants. 

For  this  problem,  the  transfer  funcdons  are  parddoned  as 

and  the  controller  is  expressed  as 

ti(5)  =  Gc(s)y(s)  (4.25) 

The  closed-loop  transfer  function  is  then  given  by 

Tzw(s)  =  G.„(s)  +  G,„{s)G,(s)[I  -  Gy„(s)Ge(s)]-*Gj.w(s)  (4.26) 
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The  design  objective  is  then  to  determine  G,(s)  such  that  the  closed-loop  system  is  stable  and 
meets  the  performance  ^ecificadons.  /  Although  general  soludons  Imve  been  obtamed  for  finite 
dimensional  systems,' the  simadon  for  infinite  dimensional  systems  is  much  more  cbmpHcated.  As 
a  result^  only  simple  systems  have  been  considered.  For  example,  the  coprime  factorizadon 
technique  is  applied  to  a  single  totsional  element  in  Piche  (1986b).  The  extension  of  such  an 
approach  to  complex  structures  would  indeed  be  a  significant  achievement. 

It  is  assumed  that  the  controller  has  fimte  order,  so  that  it  can  be  physically  implementable. 
The  n-th  order  controller  has  the  state-space  form 

x,.(x,t)  =  A,.x,.(ty+  B j'ft) ,  u(t)  =  +  D^yCO  (4.27a, b) 

which  is  represente4  in  the  frequency  domian  by 

Gc(s)  =  C,(sI.A,)->B,  +  D.  (4.28) 

Tne  objecrive,  then,  is  to  find  the  matrices  A^,  B^,  Cj,  and  Dj,  that  both  stabilize  the  closed-loop 
system  and  minimize  Tjw(s)  in  some  sense.  A  method  of  selecdng  the  order  of  the  controller  is 
also  required.  The  only  data  available  are  the  partitioned  transfer  function  matrices,  which  are 
mathemarically  exact  at  all  frequencies.  The  optimal  solution  would  then  be  valid  for  the  exact 
mathematical  model,  rather  than  some  truncation  of  it.  As  a  result,  the  modelling  error  is  restricted 
to  the  delation  of  the  mathematical  model  from  the  actual  physical  structure.  This  will  result  in  a 
less  conservative  control  design  approach  and,  consequently,  enhanced  performance. 

4.3  Limitations  of  the  TEM  Control  Design  Methodology 

Although  the  control  designs  based  on  the  TEM  methodology  have  demonstrated 
remarkable  perfoimance  (at  least  in  the  open-loop  case),  it  is  important  to  note  some  limitations  of 
this  approach.  First,  the  control  actuation  is  available  only  on  the  boundaries  of  the  structural 
elements.  For  a  small  number  of  actuators,  this  may  be  overcome  by  dividing  each  element  into 
smaller  sub-elements  at  the  point  of  control  actuation.  If  many  actuators  are  employed,  however. 
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Ais  ^proach  is'clearly  not  fduible.  AnoAer  limitadbn  is  the  requirement  that  the  initial  conditions 
'  be  zero.'  Only  in  this  way  was  it  possible  to  obtain  a  simpli  expression  for  the  control  solution. 
Nonzero  imrial  conditions  introduce  ah  extra' term,  qj(s),  in'the  dynamic  srififness  equation  for  each 
strucrural  element.  In  addition  to  making  the  optimal  control  expressions  more  complex,  this  term 
must  be  computed  by  integrating  over  the  domain  of  the  element,  as  described  in  See.  2.3.1. 
Consequently,  the  treatment  of  initial  conditions  (and  disnibmed  forcing,  for  that  matter)  increases 
the  computation  time  associated  with  the  TEM  approach  consideiably,  as  a  numerical  integration  is 
required  for  each  element  at  each  complex  frequency  of  interest. 

Some  of  these  limitations  can  be  overcome  by  working  with  the  original  PDE  for  the 
element,  expressed  in  the  time-domain.  This  forms  the  basis  of  the  direct  PDE  modelling 
approach,  which  leads  naturally  to  a  different  type  of  control  theory.  The  direct  approach  is  the 
subject  of  the  next  two  chapters. 
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5  DIRECT  PDE  MODELLING 

All  physical  systems  are  distributed  in  nanire.  This  fact  is  a  consequence  of  the  laws  of 
ptiyrics,  which  always  take  the  fonn  of  a  set  of  field  equadons  which  must  be  sadsfied  over  each 
infinitestifflal  region  in  the  spadal  domain  of  interest  As  a  result  any  exact  model  of  a  physical 
system  must  be  of  infinite  order.  Lumped  parameter  models  are,  in  general,  low-frequency 
approximadons  to  these  field  equadons.  Examples  include  lumped  electrical  component  models 
(such  as  capacitors,  resistors  and  inductors),  rigid  body  structural  idealizadons,  and  finite  element 
models.  In  this  last  example,  the  finite  order  approximadon  is  achieved  by  resiricdng  the 
deformadonal  degrees  of  freedom  of  the  system  rather  than  employing  a  low-frequency 
approximadon  directly,  but  the  result  is  essenaaily  the  same:  The  model  fails  to  recover  the  high- 
frequency  dynamics  of  the  system.  In  this  chapter,  we  introduce  the  concept  of  a  distributed, 
infinite-order  model  of  a  system,  which  retains  the  dynamics  of  the  physical  system  at  all 
frequencies.  This  approach,  hereinafter  refetred  to  as  the  direct  PDE  modelling  approach,  is 
superior  to  the'TEM  approach  when  forces  of  a  distributed  nature  act  within  the  spadal  domain  of 
the  structural  elements.  Such  forces  include  aerodynamic  and  gravitadonal  loads,  inerdal  forces, 
and  distributed  control  actuators. 

Distributed  system  models  can  be  characterized  in  either  of  two  forms.  The  first  is  an 
integral  form,  in  which  the  response  of  the  system  at  a  pardcular  drae  is  determined  by  integradng 
(with  respect  to  time  and/or  space)  the  product  of  the  distributed  forcing  inputs  and  a  Green's 
funedon  kernel.  Here,  the  Green's  funedon  relates  the  response  of  the  system  at  some  arbitrary 
point  and  time  to  an  impulse  applied  at  some  other  point  and  dme.  Thus,  this  characterizadon  is 
global  in  nature.  Given  this  approach,  it  is  possible  to  develop,  for  example,  a  distributed  control 
theory.  The  work  of  Brogan  (1968)  proceeds  along  these  lines.  However,  the  Green's  funedon 
for  an  arbitrary  system  is  extremely  difficult  to  obtain.  Indeed,  analydeal  expressions  are  only 
available  for  the  simplest  of  cases.  The  other  characterizadon  is  differendal  in  nature.  Here,  partial 
differential  equations,  describing  the  local  behavior  of  the  system,  are  used  to  develop  a  system 
model.  This  characterization  is  mush  easier  to  obtain,  as  the  physical  laws  that  describe  the  system 


are  always  local.  Consequently,  more  emphasis  has  been  placed  bn  developing  the  differential 
appoach  for  control  system  design.  'Breakwell  (1981),  for  example;  uses  the  differential 
descnpdori  to  obtain  solutions  to  the  boundary  control  of  a  simple  flexible  system.  The  dperential 
deseripdon  of  distributed  systems  wili  be  used  here  throughout;- 

5.1  One-dimensional  elements 

For  a  smjcmral  system  undergoing  small  deformadons,  the  underlying  differendal 
equadons  are,  of  course,  the  equadons  of  elasdcity.  A  completely  rigorous  and  exact  linear 
structural  model  must  therefore  account  for  general  three-dimensional  defotmadon.  Hov/ever,  for 
long,  slender  strucniral  elements,  the  defotmadon  is  primarily  a  funcdon  of  posidon  along  the 
element.  The  variadon  in  defotmadon  with  respect  to  the  other  two  directions  can  usually  be 
expressed  in  terms  of  the  deformation  along  the  length  of  the  element.  Therefore,  only  one  spatial 
coordinate  is  needed  to  describe  the  dynamics.  This  is  the  basis  for  the  axial  rod  and  Bemoulli- 
Euler  beam  models  discussed  in  Chapter  2.  While  these  idealizations  fail  to  hold  at  extremely  high 
frequencies,  their  ranges  of  validity  are  much  greater  than  those  of  finite-order  representations, 
such  as  finite  element  models. 

5.1.1  General  Formulation 

We  will  restrict  our  attention  to  one-dimensional,  linear,  time-invariant  distributed  systems. 
Such  systems  can  be  written  in  the  form 

x(x,t)  =  Lx(x)x(x,t)  +  B;^(x)u(x,t)  +  Dx(x)n(x,t) ,  xs[0,l],  ts[0,~)  (5.1) 

where  x  is  the  state  vector,  u  is  the  distributed  control  input,  and  n  is  the  distributed  disturbance 
input.  In  contrast  with  lumped-parameter  state  space  models,  these  vectors  exhibit  both  spatial  and 
temporal  dependance.  Also,  L^,  B^  and  D,;  are  linear  ^possibly  spatially  varying)  matrix 
operators.  Note  that  the  spatial  domain  has  been  normalized  to  unity.  The  boundary  conditions  are 
assumed  to  be  homogeneous,  and  ate  expressed  as 


89 


x(0,t)  =  x(l,t)  =  0,  ts[0,~) 


(5.2) 


The  entire  development  presented  here  applies  to  modelling  of  systems  with  homogeneous 
boundary  conditions  only.  It  is  therefore  assumed  that  no  control  or  disturbance  forces  are  applied 
at  the  boundaries  of  the  system,  finally,  the  inirial  conditions  are  expressed  as 

x(x,0)  =  Xo(x),  xe[0,l]  (5.3) 

5.1.2  Bemoulli-Euler  Beam  Example 

One  of  the  simplest  examples  of  a  onenlimensional  distributed  parameter  system  is  a 
Bemoulli-Euler  beam.  A  diagram  of  the  physical  system  is  shown  in  Fig.  5-1.  The  requirement 
that  the  boundary  conditions  for  the  mathematical  model  be  homogeneous  corresponds  to  pinned- 
pinned  boundary  conditions  for  the  beam,  as  will  be  shown  in  the  next  subsection.  In  addition  to 
casting  the  equation  of  motion  of  the  beam  in  the  form  given  by  Eq,  (5.1),  the  following 
subsections  describe  a  method  for  simulating  the  response  of  the  beam  system  to  various  control 
and  disturbance  forces. 

5. 1.2.1  Normalization  of  Equations  of  Motion 

In  dimensional  form,  the  beam  dynamics  are  described  by 

^[EI(x)|j^V(i(xd.t<i)l  +  m(x)|TV(i(X(i,t(i)  =  fi(X(i,t(j),  X(ie[0,L]  (5.4) 

*dL  J  *d 

where  x^i  and  tj  are  the  dimensionalized  spatial  and  temporal  variables,  respectively,  vj  is  the 
transverse  deflection,  f^  is  the  applied  distributed  control  and/or  disturbance  force,  L  is  the  beam 
length,  EI(x)  is  the  bending  stiffness,  and  m(x)  is  the  beam  mass  per  unit  length.  To  this  equation 
we  must  add  the  boundary  conditions 

Vd(('*)  =  ^Vd(0-fd)  =  Vd(l,t(i)  “  ^Vd(l,ti)  =  0  (5.5) 

dXj  SXji 
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Fig,  5-1:  Bemoulli-Euler  beam  models:  (a)  Unifcrm  beam,  (b)  Tapered  beam. 


and  the  initial  conditions 


Vd(xd.O)  “  ^Vd()Cd,0)  =  vc^jCxd)  (5.6a, b) 

We  now  introduce  nondimensional  independant  variables  according  to 

=  (5.7a.b) 

and  the  normalized  deflection  and  distributed  force  as 

v(x,t)  =  ^Vd(Xd,td),  f(x,t)  =  ^fd(Xd,td)  (5.8a, b) 

We  can  also  parametrize  the  bending  stiffness  and  mass  per  unit  length  by 

where  and  P  are  nondimensional  funcdons.  These  normalizadons  lead  to  the  following 
nondimensional  form  of  the  equadon  of  modon 

p[n(x)pv(x,t)]  +  p|5|^v(x,t)  =  f(x,t)  =  f„(x,t)  +  fn(x,t)  (5.10) 


Here,  and  fn  represent  the  normalized  distributed  control  and  dismrbance  forces,  respecdvely. 
To  obtain  the  state  space  representadon  of  the  dynamics,  we  define  the  state  vector  and  control  and 
disturbance  scalars  by 


x(x,t)  = 


tl(x)— Z^Cx,!) 

QX 


U(x,t)  =  f„(x,t) ,  n(x,t)  =  fn(x,t) 


(5.11a-c) 


The  first  and  second  elements  of  x  correspond  to  the  normalized  curvature  and  velocity  of  the 
bending  modon,  respectively.  These  choices  for  the  state  vector  components  ensure  the  well- 
posedness  of  the  system  model,  as  explained  by  Richtmyer  (1957).  Equadon  (5.5),  which  is  a 
consequence  of  the  pinned-pinned  boundary  conditions,  ensures  that  x(0,t)=x(l,t)=0  for  all  values 
of  L  The  equadon  of  modon  then  takes  the  form  of  Eq.  (5.1),  with 
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LxW  =  •’itW  =  dx(x)  = 


(5.12a.b) 


5.1.2.2  The  Case  of  Curvature  Actuation 

In  most  active  structural  control  applications,  it  is  difficult  to  implement  lightweight  inertial 
force  actuators.  This  is  particularly  challenging  for  space-based  structures,  where  stringent 
constraints  are  placed  on  strucnrral  mass.  As  a  result,  practical  structural  control  actuators  are 
usually  imbedded  within  the  structure  itself,  and  are  capable  of  producing  only  relative  deformation 
between  points  on  the  structure.  For  example,  deLuis  (1989)  demonstrates  how  an  embedded 
piezoelectric  acmator  can  be  used  to  induce  a  local  curvature  in  the  beam.  Many  such  acmators, 
placed  along  the  span  of  a  large,  beam-like  structure  wDl  then  approximate  distributed  curvature 
actuation. 

It  is  therefore  useful  to  develop  the  model  of  a  beam  with  a  distributed  curvamre  acmator. 
Such  is  the  limiting  case  of  a  beam  with  many  embedded  piezoelectric  actuators  distributed  along 
its  span.  For  this  system,  the  equation  of  motion  is  modified  to 

§[Tl(x)^v(x,t)]  +  p^gv(x,t)  =  |^m„(x,t)  (5.13) 

where  mu(x,t)  represents  the  net  action  of  the  distributed  piezoelectrics.  The  state  vector,  x,  is 
unchanged,  as  is  but  u  and  b;,  must  be  modified  to 

u(x,t)  =  mu(x,t).  b;^  =  [p(,)]0  (5.14a,b) 


5. 1.2.3  Numerical  Simulation  Using  Laplace  Transform 

Given  the  beam  dynamics  model,  there  remains  the  problem  of  actually  simulating  the 
response  of  the  beam  to  control  and  disturbance  forces.  Various  methodologies  exist  to  achieve 
this  end.  At  one  extreme,  the  dynamics  equation  is  discretized  in  both  space  and  time  and  then 
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integrated  forward  in  tinie.  This  constitutes  a  partial'differential  equation  with  ntixed  boundary  and 
Lniti^  conditions.  Although  this  method  is  widely  used,  it  requires  rather  fme  discretizations  in 
both  the  temporal  and  spatial  dimensions  to  achieve  accurate  results,  and  errors  tend  to  accumulate 
in  time.  At  the  other  extreme,  one  can  Laplace-transform  the  equation  into  the  s-domain  and  search 
for  analytic  solutions.  However,  due  to  the  distributed  nature  of  the  control  and/or  disturbance 
forces,  this  transformation  results  in  an  integro-partial  differential  equation  rather  than  a  simple 
ordinary  differential  equation  (as  would  be  the  case  for  boundary  forcing  only).  Due  to  the 
generality  of  the  distributed  forces,  a  general  analytical  solution  is  not  available. 

In  order  to  achieve  accurate  solutions  with  relatively  coarse  disaetizations,  a  third 
alternative  is  proposed.  The  dynamics  equation  is  Laplace-transformed,  resulting  in  the  above 
mentioned  integro-partial  differential  equation.  At  each  desired  complex  frequency,  a  finite 
differencing  scheme  is  used  to  solve  for  the  displacement  field.  The  data  from  a  set  of  frequencies 
is  collected,  and  the  numerically  robust  inverse  Laplace  transform  algorithm  described  in  Section 
2.2  is  used  to  conven  the  data  back  into  the  time-domain.  Because  the  transformed  equation 
represents  a  boundaiy  value  problem,  it  is  anticipated  that  its  approximate  solution  will  be  more 
stable  and  accurate  than  the  corresponding  solution  to  the  mixed  problem  associated  with  time- 
domain  integration.  The  stability  and  accuracy  of  the  inverse  transform  algorithm  has  already  been 
demonstrated  in  Chapters  2  and  3. 

The  development  presented  here  corresponds  to  distributed  force  actuation  only,  and  the 
case  of  distributed  curvature  actuation  is  addressed  in  Appendix  B.  We  first  transform  Eq.  (5.10) 
into  the  frequency-domain: 

“  fu(x,s)  +  fn(x,s)  (5.15) 

The  normalized  frequency,  s,  is  related  to  the  dimensional  frequency,  sj,  by 

s  =  SdV?”  (5.16) 
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In  this  way,  we  can  relate  the  transform  pair  v(x,t(i)  ~  vfx.s^)  with  the  pair  v(x,t)  ~  v(x,s).  We 
now  assume  that  the  feedback  control  law  is  distributed  and  linear,  and  relates  the  components  of 
the  state  vector  to  the  control  input  by 
1 

fu(x.O  =  -J  k(x,y)'^x(y,t)  dy 
1 

=  -  J  [ki(x,y)Ti(y)^v(y,t)  +k2(x,y)  |^v(y,t)j  dy  (5.17) 

0 

In  this  last  equation,  ki(x.y)  and  k2(x,y)  represent  feedback  gain  kernels.  This  type  of  control  law 
will  arise  in  the  next  chapter  as  the  optimal  solution  to  the  distributed  control  problem.  Substituting 
the  Laplace-transformed  version  of  this  control  law  in  Eq.  (5.15)  leads  to 

1 

v(x,s)]  +  ^v(x,s)  +  J  [ki(.x,y)ii(y)|3v(y,s)  +  s  k2(x,y)  v(y,s)j  dy 
0 

=  fn(.V)  +  [vo(x)  +  svo(x)]  +  Jk2(x,y)  vo(y)  dy  (5.18) 

The  term  involving  kj  in  this  equation  can  be  integrated  by  pans  twice  so  that  the  derivative  with 
respect  to  y  operates  on  kj.  The  boundary  term  arising  from  this  operation  vanishes  due  to  the 
homogeneous  boundary  conditions.  By  making  the  following  associations 


k(x.y,s)  =  ^[ki(x,y)ii(x)]-i-sk2(x,y) 

(5.19a) 

f(x,s)  =  fn(x,s)  +  fi(x,s)  +  f<;(X.S) 

(5.19b) 

fi(x,s)  =  ^  [vo(x)  +  s  vo(x)] 

(5.19c) 

1 

fc(x,s)  =  1  k2(x,y)  vo(y)  dy 

(5.19d) 

0 


the  dynamics  equation  reduces  to 
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(5.20) 


1 

^[ti(x)^v(x,s)]  +  ^v(x,s)  +  J  k(x,y,s)  v{y,s)  dy  =  f(x,s) 

A  similar  result  is  avaDable  for  the  case  of  curvature  actuadon,  and  can  be  found  in  Appendix  B. 

Equation  (5.20)  must  be  solved  numerically  for  each  value  of  s  needed  to  construct  the  time 
response.  To  do  so  requires  a  discretization  of  the  spatial  domain  into  N  uniform  subregions.  The 
boundaries  of  these  subregions  are  given  by 

Xi  =  ^.  i  =  0 . N  (5.21) 

We  can  now  use  the  values  of  f(x,s)  evaluated  at  these  Xj  to  determine  v(x,s)  at  these  same 
coordinates.  By  defining  the  vectors 

v(s)  =  {v(xi,s))  .  f(s)  »  {r(Xi.s)}  (5.22a.b) 

and  the  matrix 

K(s)  »  [k(xi,yj.s)]  (5.23) 

an  approximation  to  Eq.  (5.20)  is  easily  obtained.  The  first  term  is  replaced  by  the  finite  difference 
approximation 

0[il(x) § v(x.s)]  »  DHD  v(s)  (5.24) 


where  D  is  a  constant  banded  matrix  of  coefficients  representing  the  second  derivative  operation: 


r2-2 

-1  2-1  0 


0  -1  2-1 
-2  2-J 


and  H  is  the  discretized  representation  of  il(x): 

H  =  diag  [iKxj)] 


(5.25) 


(5.26) 
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The  second  term  in  Eq.  (5,20)  is  trivially  apprwdmated  by 

^v(x.s)  -  s2Bv(s)  (5.27) 

where 

B  .  dl..[p^]  (5.28) 

Finally,  the  integral  term  is  replaced  with  a  summadon: 

1 

Jk(x,y,s)v(y,s)dy  »  ^K(s)v(s)  (5.29) 

Collecting  terms,  the  discrerized  equadon  becomes 

[•N4  D  H  D  +  s2  B  +  i  K(s)]  v(s)  =  f„(s)  +  f,(s)  +  fc(s)  (5.30) 

Thus,  a  single  matrix  inversion  is  required  at  each  complex  frequency.  If  the  frequencies  required 
for  the  inverse  Laplace  transform  are  given  by 

s  =  Si . Sn  (5.31) 

then  the  solutions  of  Eq.  (5.30)  can  be  grouped  according  to 

V  =  {;(si)  ...  v(s„)}  (5.32) 

The  time-domain  responses  at  each  Xj  arc  then  obtained  by  applying  the  inverse  transform 
algorithm  to  each  row  of  V. 

Figure  5-2  presents  the  response  of  a  uniform  and  a  linearly  tapered  beam  to  a  sinusoidal 
ininal  displacement  and  zero  initial  velocity.  (The  plots  display  time  and  x-coordinate  along  the 
beam  as  independent  variables,  with  transverse  deflection  along  the  vertical  axis.)  For  the  uniform 
beam,  these  initial  conditions  correspond  to  the  second  mode  of  vibration.  Consequently,  no  other 
mode.s  are  excited,  as  can  be  seen  in  the  figure.  For  the  tapered  beam,  other  modes  become 
involved,  as  the  individual  mode  shapes  are  more  complicated.  Hgure  5-3  displays  the  simulation 
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beam,  (b)  Tapered  beam  (Beam  diameler  varies 


results  for  a  unifomi  beam  impacted  with  a  unit  tiansveise  impulse  at  its  center-span.  The  plot  on 
the  left  corresponds  to  a  long  time  scale,  and  indicates  that  the  resulting  modon  is  predominantly 
composed  of  a  first  mode  vibradon  with  odd  harmonics.  The  plot  on  the  right,  which  corresponds 
to  a  shorter  time  sede,  accentuates  the  wave-like  characterisdes  of  the  response.  The  disturbance, 
which  begins  at  the  center-span,  quickly  moves  towards  the  boimdaries  and  reflects.  These 
reflecdons  evenmally  set  up  the  complex  modal  modon  observed  in  the  long  dme  scale  plot.  Note 
that,  for  a  short  dme  following  the  impact,  the  deflecdon  of  the  center-span  varies  as  the  square- 
root  of  time.  This  behavior  agrees  with  the  beam  theory  presented  by  Nowacki  (1963),  where  the 
response  of  a  beam  of  infinite  extent  is  addressed.  Note  also  that  the  disnirbance  reaches  the 
boundaries  almost  instandy,  which  is  characterisdc  of  the  Bemoulli-Euler  beam  assumprion  of  no 
cross  secdon  rotary  inertia.  This  effect  is  more  apparent  in  Fig.  5-4,  where  the  responses  of  a 
Bemoulli-Euler  beam  (on  the  left)  and  a  Timoshemko  beam  (on  the  tight)  are  compared.  The 
simuladon  of  the  Timoshenko  beam  is  discussed  in  Appendix  C.  For  these  simuladons,  free-free 
boundary  condidons  are  assumed,  and  the  impact  occurs  at  a  boundary.  The  effect  of  rotary  inteda 
is  immediately  apparent,  and  manifests  itself  as  a  finite  disnirbance  propagadon  velocity  in  the 
beam.  Also,  the  reflection  of  the  shear  wave  p.'opagating  through  the  Timoshenko  beam  can  be 
seen  in  the  plot  on  the  right  of  the  figure. 

5.2  Two-dimensional  elements 

Many  element  models  require  two  independant  spadal  coordinates  to  specify  the  domain  of 
the  element.  These  models  include  membranes,  plates  in  bending,  shells,  and  plane  stress 
elements.  In  all  cases,  the  third  spatial  dimension  is  of  sufficiently  small  extent  in  comparison  with 
the  other  two  dimensions  so  that  a  two-dimensional  idealization  is  reasonably  accurate.  These 
models  have  the  general  differential  form 

x(x,y,t)  =  Lxx(x,y,t)  +  BxU(x,y,t) ,  x,ye[0,l],  ts[0,~>)  (5.33) 

Two  examples  of  two-dimensional  elements  are  given  below. 
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5.2.1 .  Membrane  model 

The  normalized  equation  of  motion  of  a  membrane  is  given  by 

-V2v(x,y.t)  +  v(x.y,t)  =  f(x,y,t)  (5.34) 

where  v  represents  the  normalized  deflection,  f  represents  a  normalized  force  per  unit  area,  and 


! 

is  the  Laplacian  operator.  Taldng  the  deflection  and  its  velocity  as  the  state  variables: 

1 

x(x,y,t)  =  .  .  u(x,y,t)  =  f(x.y,t) 

Lv(x,y,t)J 

(5.35a.b) 

1 

and  defining  the  operators  and  b;;  by 

1 

(5.36a, b) 

1 , 

1 

leads  to  the  relation  given  by  Eq.  (5.33). 

r 

5.2.2  Plate  model 

Another  two  dimensional  element  is  a  plate  in  bending.  Here,  the  equation  of  motion  is 

i. 

V‘*v(x,y,t)  +  v(x,y,t)  =  f(x,y,t) 

(5.37) 

i; 

In  this  case,  the  following  state  vector  and  forcing  input  definitions  are  appropriate; 

[ 

r7^v{x,y,t)l 

x(x,y,t)  =  .  ,  u{x,y,t)  =  f(x,y.t) 

L  v(x,y,t)  J 

(5.38a, b) 

i 

The  equation  of  motion  then  takes  the  form  of  Eq.  (5.33),  with 

l 

1 

w-[vn.  ■>.-[;] 

(5.39a,b) 

1 

1 
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5.2.3  Complexity  Issues 

Unfortunately,  the  direct  simulation  of  two-dimensional  elements  is  considerably  more 
computadonally  intensive  than  the  simulation  of  one-dimension  elements.  Because  the  spadal 
domain  is  given  by  two  independant  variables,  discredaadon  in  two  dimensions  in  required.  As  a 
result,  no  simuladon  results  are  cuirendy  available.  A  complete  development  of  the  direct 
siffiuladon  of  two  dimensional  elements  is  the  subject  of  future  research. 

5.3  Multiple  Element  Formulation 

Distributed  modelling  and  simuladon  of  muldple  member  strucnires,  such  as  space 
hrames  and  trusses,  is  a  considerably  more  difficult  problem  than  the  single  element 
siniarion,  even  for  one-dimensional  elements.  The  primary  difficulty  is  in  the  mathemadcal 
treatment  of  the  boundary  conditions  that  arise  at  element  junctions.  A  rigorous,  general 
assembly  procedure  for  complex  structures  using  direct  PDE  modelling  remains  to  be 
developed.  One  approach  currently  considered  is  to  define  a  normalized  local  coordinate 
system  (x=0  at  one  end  of  an  element  and  x=l  at  the  other  end)  for  each  structural  member, 
as  shown  in  Fig.  5-5,  and  collect  the  states  associated  with  each  element  into  a  large  state 
vector.  The  dynamics  of  the  entire  structure  is  then  still  represented  by  Eq.  (5.1),  and 
becomes  block  diagonal,  with  each  block  representing  the  dynamics  of  one  member.  Tlie 
boundaty  conditions  then  relate  various  elements  of  the  state  vector  at  x=0  and/or  x=l.  The 
difficulty  in  the  direct  modelling  approach  then  lies  in  utilizing  these  awkward  boundary 
conditions  for  the  purposes  of  simuladon  and  control  design.  The  direct  multiple  clement 
formulation  remains  an  open  area  of  research. 
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6  CONTROL  DESIGN  BASED  ON  DIRECT  PDE  MODELS 

This  chapter  studies  the  properties  of  optimal  sol'  idons  to  the  distributed  control  of  systems 
described  by  direct  PDE  models.  By  distributed  control,  we  imply  that  control  effort  is  imparted  to 
the  structure  in  a  spatially  distributed  sense.  This  can  be  thought  of  as  the  limiting  process  of 
employing  an  increasing  number  of  actuators,  all  of  which  have  a  decreasing  spatial  domain  of 
influence.  Particular  emphasis  is  placed  on  a  specific  example,  that  of  a  Bemoulli-Euler  beam. 
Both  fuilte-  and  infinite-length  beam  systems  are  snidied,  and  comparisons  are  made  between  the 
corresponding  optimal  solutions.  Also,  two  types  of  beam  actuation  are  addressed.  The  first  is 
force  actuation,  which  is  commonly  used  in  theoretical  studies  yet  is  rarely  achievable.  The  second 
is  curvature  actuation,  which  is  more  realizable  (as  mentioned  in  the  previous  chapter)  but  less 
often  addressed  in  theoretical  works. 

Distributed  control  is  by  no  means  a  new  topic.  In  fact,  the  essential  mathematical 
groundwork  was  established  in  the  1960's  by  Butkovskii  (1960),  Wang  (1964),  and  Lions 
(1971).  The  results  then  obtained  were  analogous  to  the  classical  LQR  solution  (e.g.,  the  Riccati 
matrix  equation  was  replaced  by  a  Riccati  operator  equation),  and  were  derived  using  the  principle 
of  optimality  and/or  advanced  functional  analysis.  A  later  work  by  Tzafestas  (1970)  derived  the 
necessary  conditions  for  optimatity  from  a  variational  calculus  approach.  A  mathematically 
rigorous  derivation  of  the  Riccati  operator  equation  is  performed  by  Gibson  (1979).  Also,  Balas 
(1982)  addresses  several  implementation  issues,  including  the  use  of  a  finite  set  of  sensors  and 
acniators  and  a  finite  order  controller.  Perturbation  methods  are  utilized  to  determine  criteria  for 
closed-loop  stability.  Until  now,  the  complexity  of  the  distributed  control  problem  has  rendered  it 
a  mere  mathematical  curiousity,  rather  than  a  practical  tool.  With  today's  computer  resources, 
however,  the  solutions  to  simple  problems,  such  as  the  Bemoulli-Euler  beam  system  described 
below,  are  within  reach. 


105 


6.1.  Linear  Quadratic  Optimal  Control  Theory,  in  the  1-D  Case 

The  Lheoiy  presented  in  this  subsection  closely  resembles  the  work  of  Tzafestas  (1970). 
We  restrict  our  attention  to  one-dimensional,  lin^,  time-invariant  distributed  systems.  Such 
systems  can  be  written  in  the  form 

x(x,t)  =  Lx(x)x(x,t)  +  Bx(x)u(x,t) ,  xe  [0,1] ,  t€  [0,~)  (6.1) 

which  is  identical  to  Eq.  (5.1),  except  that  the  disnirbance  input  has  been  set  to  zero.  As  before, 
the  boundary  conditions  are  assumed  to  be  homogeneous,  and  are  expressed  as 

x(0,t)  =  x(l,t)  =  0,  ts[0,~)  (6.2) 

while  the  initial  conditions  are  expressed  as 

x(x,0)  =  xo(x),  xs[0,l]  (6.3) 

The  simplest  optimal  distributed  compensator  is  derived  under  the  assumption  of  full  state 
feedback.  That  is,  perfect  measurements  of  x(x,t)  are  available  in  a  continuous  sense  throughout 
the  spatial  domain  at  eveiy  instant  of  time.  It  is  also  assumed  that  control  actuation  is  available  in  a 
similar  distributed  sense.  While  these  assumptions  are  rather  crude,  they  serve  to  define  an  upper 
limit  of  achievable  performance  for  the  control  system  to  be  designed.  The  optimal  distributed 
control  problem  can  then  be  stated  as  follows:  Given  an  arbitrary  initial  condition,  detemune  the 
control  requhed  to  return  the  system  to  the  zero  state  while  minimizing  some  cost  criterion.  We 
will  assume  a  linear  quadratic  cost  functional  of  the  form 

“1 

I  =  jJJ  [x(x,t)TQ(x)x(x,t)  u(x,t)TR(x)u(x,t)]  dx  dt  (6.4) 

where  Q  and  R  are  symmetric  (possibly  spatially  varying)  weighting  matrices.  This  is  the 
distributed  analogue  of  the  classical  linear-quadratic  regulator  (LQR)  problem.  Its  solution  is 
obtained  by  extending  the  classical  variational  calculus  approach  to  distributed  systems.  Note  that 


the  cost  functional  has  infinite  dote  hoiizon.  This  corresponds  to  the  steady-state  LQR  problem. 
The  finite  time  problem  is  also  of  interest,  but  provides  no  addidonal  insight 

We  first  augment  the  cost  funcdonal  with  the  system  dynamics  via  a  costate  vector,  p(x,t): 

ool 

Ja  =  J  +  JJ  p(x,t)‘^[Lx(x)x(x,t)  +  Bx(x)u(x,t)  -  x(x,t)]  dx  dt  (6.5) 

The  augmented  cost  funcdonal  now  depends  on  the  three  vectors,  x,  u,  and  p.  The  cost  is 
minimized  by  setting  the  variation  in  cost  due  to  independent  perturbations  in  these  three  vectors  to 
zero.  Thus, 


■>1 


5la(Sp)  =  ii5p(x.t)‘^[  Lj((x)x(x,t)  T  B;j(x)u(x,t)  -  x(x.t)]  dx  dt  =  0 ,  V  5p(x,t)  (6.6a) 


00 


~1 

5Ja(Su)  =  j|  [u(x,t)^R(x)5u(x.t)  +  p(x,t)'*^Bx(.x)5u(x,t)]  dx  dt 
00 


JJ  [u(x,i)^R(x)  +  (Bx(x)p(.x,t)]''^]5u(x,t)dxdt  =  0,  V  5u(x,t)  (6.6b) 


eol 

5Ja(5.x)  =  [ x(x,t)TQ(x)Sx(x,t)  +  p(x,t)'’^Lx(x)6x  -  p(x,t)^6.x(x,t)]  dx  dt 


©ol 

=  JJ  [x(x,t)TQ(x)  +  [Lx(x)p(x,t)]'^  +  p(x,t)'^]6x(x,t)dxdt  =  0,  V  5x(x,t)  (6.6c) 
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Equation  (6.6a)  recovers  the  system  dynamics,  while  Eq.  (6.6b)  determines  the  control  law.  The 
superscript  (*)  represents  the  formal  adjoint  operator,  which  is  defined  by 


1  1 

J  a(xfL;,(xMx)  dx  =  j  [L;;(x)a(x)]Tb(x)  dx  (6.7) 

for  homogeneous  boundary  condidons.  For  a  linear  spadal  differendal  operator,  its  adjoint  is 
detenruned  by  integrating  by  parts  with  respect  to  the  spatial  dimension.  Solving  for  u  in  Eq. 
(6.6b)  yields 

u(x,t)  =  -  R(x)‘’B*(x)p(x,t)  (6.8) 

The  integrated  terms  resulting  from  the  integration  by  parts  in  Eq.  (6.6c)  vanish,  due  to  the 
homogeneous  boundary  conditions  and  the  added  requirement  that 

P(0,t)  =  p(l,t)  =  0,  ts[0,«)  (6.9) 

The  third  term  in  the  integrand  of  Eq.  (6.6c)  is  integrated  by  pans  with  respect  to  time.  The 
integrated  terms  go  to  zero  due  to  the  specification  of  the  initial  conditions  for  the  system  and  the 
requirement  that 

p(x,~)  =  0,  X6(0,l]  (6.10) 


Equations  (6.1),  (6.8),  and  the  integrand  in  Eq.  (6.6c)  lead  to  the  following  equations: 

x(x,t)  =  Lx(x)x(x,t)- B;((x)R(x)'*B’(x)p(x.t) 

P(x,t)  =  -  Q(x)x(x,t)  -  Lx(x)p(x.t) 


Equations  (6.1  la,b)  represent  the  state-costate  equations  for  the  distributed  control 
problem.  Lions  (1971)  shows  that  there  exists  a  relation  between  the  state  and  the  costate  of  the 
form 

p(x,t)  =  P^(x)x(x,t)  (6.12) 
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where  is  some  linear  matrix  operator  on  x.  Substituting  this  form  in  Eq.  (6.11)  results  in  a 
nonlinear  matrix  Riccad  operator  equation  in  P^.  Such  an  equation  is,  in  general,  difficult  to  solve. 
Several  approximate  solution  techniques  are  described  in  Juang  (1983),  Schaechter  (1982),  and 
2;ambettalds  (1989).  However,  it  is  possible  to  express  the  linear  operator  in  a  different  form,  so 
that  a  solution  is  easily  attained  by  numerical  methods.  The  assumed  form  of  the  solution  is  the 
same  as  used  by  Wang  (1964)  and  Tzafestas  (1970): 

1 

p(x,t)  =  /S(x,y)x(y,t)dy  (6.13) 

0 

where  S  is  tiie  distributed-parameter  analogue  of  the  Riccati  matrix  for  lumped-parameter  systems. 
Equation  (6.9)  automatically  imposes  the  constraints 

S(0,y)  =  S(l,y)  =  0 ,  ye  [0,1]  (6.14) 

For  complete  generality,  S  must  include  generalized  functions,  such  as  Dirac  delta  functions  and 
their  derivatives,  if  necessary.  Also,  Wang  (1964)  shows  that  S  is  symmetric  in  its  arg  iments 
(i.e.,  S(x,y)  =  S(y,x)).  Using  Eqs.  (6.8)  and  (6.13),  the  feedback  control  law  becomes 
1 

u(x,t)  =  -JK(x.y)x(y,t)dy,  K(x,y)  =  R(x)-'B‘(x)S(x,y)  (6.153,b) 

Thus,  the  control  law  is  linear  and  distributed. 

It  remains  to  derive  a  relation  that  enables  the  computation  of  S.  Differentiating  Eq.  (6.13) 
and  introducing  Eq.  (6.1  la)  leads  to 
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1 

p(x,t)  =  Js(x,y)  [Ly(y)x(y,t)-By(y)R(y)'*Bj(y)p(y.t)]dy 

=  J  S(x,y)Lj(y)'''x(y,t)-S(x,y)By(y)R(y)’*Bj(y)jS(y,i)x(r,t)dzjdy 
0 

1  ^ 

=  S(x,y)L;(y)^- Js(x.z)Bj(z)R(z)-'B;(2)S(z,y)dz 

J  L  0  . 


x(y.t)dy  (6.16) 


Once  again,  an  Integration  by  parts  applied  to  the  first  term  in  the  integral  results  m  the  adjoint 
operator.  The  boundary  terms  again  vanish,  subject  to  the  restricrion 

S(x,0)  =  S(x,l)  =  0 ,  xe[0,l]  (6.17) 

It  should  be  noted  that  the  transposes  of  adjoint  operators  operate  to  the  left  in  this  case.  Similarly, 
substituting  Eq.  (6. 13)  into  the  right  side  of  Eq.  (6.1  lb)  yields 

1 

p(x,t)  =  -Q(x)x(x,t)-J  [Lj(x)S(x,y)x(y,t)]dy 
1 

=  ■  J  [Q(x)5(x.y)  +  Lj(x)S(x,y)]  x(y,t)  dy  (6.18) 

where  5(.x)  is  the  Dirac  delta  function.  Note  that  the  state  vector,  x,  is  isolated  from  each  term 
under  the  integral  in  Eqs.  (6.16)  and  (6.18).  Thus,  subtracting  these  two  equations  and  setting  the 
resulting  integrand  to  zero  yields  the  desired  relation: 

1 

L*(x)S(x,y)  +  S(x,y)  L;(y)’'  +  Q(x)  5(x-y)  -  J S(x,z)B^(z)R(z)-'Bj(z)S(z,y)  dz  =  0  (6.19) 

This  relation  is  a  functional  nonlinear  matrix  integro-partial  differential  equation  in  x  and  y,  and 
represents  the  distributed  parameter  analogue  of  the  control  algebraic  Riccati  equation.  Note  that 
we  have  assumed  S  to  be  time-invariant,  which  corresponds  to  the  steady-state  linear  quadratic 
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regulator.  For  a  finite  time  problem,  the  zero  on  the  right  hand  side  of  Eq.  (6. 19)  would  be 
replaced  by  ■§iS(x,y). 

6.2  Distributed  Control  of  a  Finite  Beam 

In  this  section,  we  apply  the  distributed  control  theory  just  presented  to  a  BemouUi-Euler 
beam  of  finite  length,  as  described  in  Lupi  (1991a).  Force  actuation  will  be  assumed  initially,  and 
curvature  actuation  will  be  deferred  to  Section  6.2.5.  The  feedback  gains  will  be  determined  by 
numerical  solution  of  the  Riccati  equations  given  by  Eq.  (6.19).  Although  these  equations  are 
quite  complex,  their  solutions  are  readily  attainable  with  the  proper  mix  of  algebraic  manipulation 
and  numerical  computation.  Pinned-pinned  boundary  conditions  are  assumed  for  the  example 
applications  presented  in  this  section. 

6.2.1  Cost  Functional 

The  dimensionalized  form  of  the  cost  functional  for  this  system  is  expressed  by 

Jd  |quWEI(x)j^|^j‘+qT(x)m(x)j^|^j^+r(x)— fj|d.X(idt(j  (6.20) 


Thus,  qu  represents  a  weighting  on  defonmational  potential  energy,  q-f  represents  a  weighting  on 
kinetic  energy,  and  r  weighs  control  effort.  The  physical  parameters  El,  m  and  L  are  introduced  so 
that  all  three  weights  have  the  same  units.  The  cost  functional  can  then  be  normalized,  yielding 


dx  dt 


(6.21) 


where  the  nondimensional  cost  is  defined  by 
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1 


(6.22) 


jr  = 


1  j 

VH(p)m(0)L2  “ 


The  cost  functional  then  takes  the  form  of  Eq.  (6.4),  with 


Q(x)  = 


qu(x)n(x)  0  ■ 
0 


.  KM  -  g 


(6.23) 


6.2.2  Derivation  of  the  Necessary  Conditions 

We  must  first  determine  the  adjoint  operators  corresponding  to  and  bjj.  This  is 
accomplished  by  using  the  formal  definition  expressed  in  Eq.  (6.7)  and  integrating  by  pans, 
yielding 

L^(.x)  =  lJ(x)  =  b*(.x)  =  bj(.x)  =  [0P(x)]  (6.24a.b) 

Using  these  expressions  in  (6.19)  yields 

1 

L:S(x,y)  +  S(x,y)L;^+Q(x)5(x-y)  -J’l^^S(x,a)  [J  ?]s(2,y)dz  =  0  (6.25) 

Also,  making  use  of  Eq.  (6. 1 5),  the  feedback  law  becomes 

«(x.t)  =  -Jk(x,y)Tx(y.t)dy,  k(x,y)  =  (6.26a,b) 

The  effect  of  the  curvature  feedback  term  (ki)  is  to  stiffen  the  beam,  which  reduces  the  settling  time 
of  the  system,  while  the  effect  of  the  velocity  feedback  term  (k2)  is  to  increase  the  damping  of  the 
system. 

Equation  (6.25)  represents  a  system  of  four  coupled,  nonlinear,  integro-partial  differential 
equations.  Due  to  the  symmetry  of  S,  the  fourth  is  redundant.  Also,  only  S  ]2  and  S22  are  needed 
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to  compute  the  feedback  gains.  The  equation  for  S12,  which  represents  the  curvanire  feedback 
gain  kernel,  is  uncoupled  &om  the  others: 

^  [PW  Si2(x,y)]  [P(y)  Si2(x,y)] 

1 

=  qu(x)tl(x)5(x-y)-J’li!E^Si2(x,z)Si2(z,y)dz  (6.27) 

Similrly,  for  the  velocity  feedback  gain  kernel,  the  relevant  integro-partial  differential  equation  is 

^2  Si2(x.y)]  +  n  [n(y)  Si2(x.y)] 

OX  oy 

1 

+  ^S(x-y)  -  J’l^^S22(x,z)S22(z,y)dz  =  0  (6.28) 

Note  that  this  second  equation  requires  knowledge  of  Si2(x,y),  which  is  determined  upon  solving 
(6.27).  Thus,  the  two  equations  must  be  solved  consecutively,  using  approximate  numerical 
methods. 

6.2,3  Numerical  Solution  of  the  Riccati  Equations 

Previous  attempts  to  obtain  a  numerical  solution  to  the  optimal  distributed  control  problem 
for  a  particular  system  have  most  often  dealt  with  the  operator  form  of  the  PJccati  equation,  which 
is  derived  by  Gibson  (1979)  using  Eq.  (6.12)  rather  than  Eq.  (6.13).  Usually,  the  solution  is 
expressed  as  a  series  expansion  of  spatial  differential  operators  of  increasing  order,  as  in  Juang 
(1983).  In  some  cases,  the  distributed  control  law  is  only  solved  at  points  where  discrete  controls 
are  to  be  applied,  which  leads  to  a  slightly  suboptimal  design.  Balas  (1982)  takes  this  appro.ach. 
However,  in  this  formulation,  the  functional  form  of  the  Riccati  equations  leads  naturally  to  a 
numerical  solution  procedure.  Because  of  the  fundamental  differences  in  the  forms  of  Eqs.  (6.27) 
and  (6.28),  a  separate  algorithm  is  developed  for  each  equation,  as  discussed  below. 
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6.2.3.1  Solution  of  the  First  Riccad  Equadoa 

Equadon  (6.27)  is  solved  by  spatially  discretizing  the  domain  of  S12  and  using  finite 
differencing  and  summation  to  approximate  the  derivative  and  integral  operations,  respectively.  A 
modified  relaxation  algorithm  is  then  invoked  to  converge  upon  the  solution.  We  begin  by 
discretizing  the  spatial  variables  according  to 

Xi  =  ^.  i=0.....N  (6.29) 

and  defining  the  mesh 

Sjj  =  Si2(.Xi.yj)  (6.30) 

A  simple  approximation  to  the  derivative  terms  is  then 

0[P(x)S,2(x.y)]  +^[P(y)S,2(x.y)]  -  [a?.  -  2(Pi+pj)Sij]  (6.31) 

where  aJ.  is  defined  by 

=  Pi-1  Si-l  j  +  Pi+l  Sitl  j  +  Pj-l  Si  j.I  +  Pj+I  Si;M  (6-32) 

The  forcing  term  in  Eq.  (6.27)  can  be  appro.ximated  by 

qu(x)Tl(x)5(x-y)  »  NqujTliSjj  (6.33) 


and  Sjj  is  the  ^screte  Kroneker  delta  function.  Finally,  the  integral  term  is  replaced  with  a 
summation,  leading  to 


J^'^^^^^Si2(x,z)Si2(z,y)dz  - 
0  k=0 


where 
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(6.35) 


N 


nkPi 


k-Q 

k»i 

k-j 


Sik  ®kj 


Note  that,  in  Eqs.  (6.31)  and  (6.34),  the  terms  involving  Sjj  have  been  isolated  from  terms 
involving  neighboring  points.  Collecting  the  approximate  expressions,  we  have,  for  the  finite 
difference  equation 

N2ri^. - 2N2(P;+Pj)Sjj  =  Nqu.Tli5jj  -  jjly  -  (6.36) 

Thus,  given  an  initial  estimate  for  the  solution  at  each  mesh  point,  s?,  the  entire  mesh  is 
successively  iterated  according  to  the  rule 


s”?’  =  S"-CDe?.,  0<Ci><2 

U  ij  u’ 


(6.37) 


In  this  last  equation,  ejl  represents  the  residual  error  at  each  mesh  point  at  the  n-th  iteration.  An 
expression  for  this  error  is  obtained  by  solving  Eq.  (6.36)  for  sy,  which  yields 

N2Af,  -k  '-Ijj  .  Nqu.qjSjj 


2N2(pi+Pj)  -  p^Sii  +  -  |5,j) 


(6.38) 


Also,  m  is  a  relaxation  parameter,  and  can  be  adjusted  to  maximize  the  rate  of  convergence  towards 
a  solution,  as  discussed  in  Press  (1986). 

The  condition  for  a  converged  solution  is  given  by 

|ei)l<e.  Vi,j  (6.39) 

where  e  is  some  small  positive  constant.  The  relaxation  method  is  guaranteed  to  converge  when  r 
approaches  infinity  (In  this  case,  Eq.  (6.27)  reduces  to  Poisson’s  equation),  and  tests  have  shown 
that  convergence  is  maintained  over  a  wide  range  of  values  of  qy  and  r,  provided  co  is  adjusted 
accordingly. 
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In  Hg.  6-1.  the  gain  surfaces  for  various  are  shown  for  the  case  of  constant  section 
properties.  The  effect  of  the  boundaries  can  be  seen  by  observing  the  gain  surface  near 
(x,y)=(0.0)  and  (x,y)=(l.l).  Qualitatively,  the  influence  of  the  boundary  conditions  extends  over 
a  smailer  domain  as  the  control  authority  is  increased  (i.e..  as  the  quantity  r/q„  becomes  smaller). 

Aquandtadve  analysis  indicates  that  the  extent  of  this  “boundary  layer"  is^^^^^ 

(r/qu)W  This  makes  sense  physically.  As  the  control  authority  is  increased  (r  decreasing),  the 
system  is  able  to  suppress  the  majority  of  a  disturbance  before  the  energy  reaches  and  reflects  off 

the  boundary.  Thus,  near  the  center  of  the  beam,  the  controller  models  the  beam  as  if  it  were 
infinite  in  length. 

6.2.3.2  Solution  of  the  Second  Riccad  Equation 

Equation  (6.28)  does  not  have  a  well-behaved  solution,  since  it  requires  that  the  integral  of 
S22(x.y)  cancel  the  delta  function.  We  therefore  ma.ke  the  foUowing  substitution: 

S22(x,y)  =  S22(x.y)  +  g(x)  5(x.y)  (6.40) 


°  I fl(x)^qu(x)  +  q^fx)]  (6.4 j) 

This  IS  equivalent  to  identifying  a  collocated  component  in  the  velocity  feedback  kernel.  Equation 
(6.28)  then  becomes 

1 

j"  S  22(x.2)S22(z,y)  dz  +  [g(x)+g(y)]  S22(x.y)  -  c(x,y)  =  0  (6.42) 

where  the  known  forcing  terra  is  given  by 

c(x,y)  =  j2[fl(x)Si2(x.y)]+0[q(y)Sj2(x,y)].q„(x)^5(x.y)  (6.43) 
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It  is  easy  to  show  that  c(x,y}  is  continuous  (assuming  ilM  PM  ^  continuous),  as  the  paitial 
differentiation  terms  produce  delta  funcdons  that  exacdy  cancel  the  term  involving  5(x-y). 

The  soludon  algorithm  for  S22  is  straightforward.  Upon  discredzing  in  the  spadal 
dimension,  Eq.  (6.42)  becomes 


$22^*^522  +  S22G  GS22  *0  =  0 

(6.44) 

§22  =  [s22(xj.yj)] 

(6.45a) 

C  =  [c(Xi,yj)] 

(6.45b) 

G  =  diag  [g(xi)] 

(6.45c) 

R  =  diagfN  — -'m 

(6.45d) 

and  N  is  the  number  of  mesh  points  between  x=0  and  x=l.  The  matrix  equation  is  solved  by 
completing  the  square.  After  ore-  and  post-multiplying  by  R'*^,  Eq.  (6.44)  can  be  factored  as 

[r-‘%2R'‘^  +  g]'  =  R-’^CR-'^  +  (6.46) 

The  right  hand  side  of  Eq.  (6.46)  is  symmetric  and  positive  semidefinite.  It  therefore  has  the 
eigenvector  decomposition 

r-1/2cr.1/2^G- = 'WAW’’'  (6.47) 

where  A  is  a  diagonal  matrix  with  non-negative  entries.  Finally,  substituting  Eq.  (6.47)  in  Eq. 
(6.46)  and  solving  for  §22  gives 

§22  =  R‘^[WA’/2vv'^-  G]r'/2  (6.48) 
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Typical  k2(x,y)  surfaces  are  shown  in  Fig.  6-2  (ic2(x,y)  is  just  k2(x,y)  without  the  delta  function 
conesponding  to  the  collocated  feedback  component}.  Finally,  Fig.  6-3  shows  the  feedback  gain 
kernels  associated  with  a  tapered  beam. 

6.2.4  The  Case  of  Curvature  Actuation 

For  the  sake  of  simplicity,  we  will  assume  constant  section  propeities  and  weighting 
functions  for  this  case.  This  makes  it  possible  to  obtain  an  analytical  expression  for  the  feedback 
gains.  G^’ote  that,  for  force  actuation,  a  numerical  solution  procedure  is  required  even  in  the  case 
of  constant  section  properties  and  weightings.) 

The  feedback  gains,  expressed  in  terms  of  the  solutions  to  the  Riccarf  equations,  become 

,  ki(x,y)  =  77^Si2(x,y),  k2(x,y)  »  7|^S22(x,y)  (6.49a, b) 

and  the  functional  Riccati  equations  themselves  become 

1 

|^Si2(x,y)  +  pSi2(x,y)  =  qu5(x-y)  •  7  JSi2(x,z)^Si2(a,y)dz  (6.50a) 

X  y  J  i 

1 

0Si2(x.y)  +  |^Si2(x,y)  +  qT5(x-y)  -  7  JS22(x,z)^S22(z-y)(3z  =  0  (6.50b) 

X  y  ^  z 

Integrating  by  pans  and  invoking  the  homogeneous  boundary  conditions  yields 

1 

0Si2(x,y)+^Si2(x,y)  =  qu5(x-y)  -  7  J  0Si2(x,z)0Si2(z,y)dz  (6.51a) 

0 

1 

0Si2(x,y)  +  pSi2(x,y)  +  qT5(x-y)  -  7  J^S22(x,z)^S22(z.y)dz  =  0  (6.51b) 

^  0 

Funhermore,  introducing  Eqs.  (6.49a, b)  and  exploiting  the  symmetry  of  S(x,y)  yields 


2rki2(x,y)  =  qoS(x-y)  -  rjki2(x,z)ki2(z,y)dz  (6.52a) 

1 

rjfc22(x.z)k22(z.y)dz  =  2rki2(x,y)  +  qi.5(x-y)  (6.52b) 

It  can  easily  be  shown  (hat  the  generalized  funcdons 

ki(x.y)  =  ["s/ 1  +  ‘’f-  -  1  ]  8(x-y)  (6.53a) 

k2(x.y)  =  V ?+  2  +  -  1  ]  5(x-y)  (6.53b) 

solve  the  Riccati  equations.  The  optimal  control  is  therefore  purely  collocated,  even  though  the 
controller  has  access  to  the  state  vector  over  the  entire  spatial  domain. 

6.2.5  Closed-Loop  Simulation  Results 

The  closed-loop  simulations  of  a  uniform  beam  with  a  sinusoidal  initial  condition  and 
various  control  and  state  weightings  are  shown  in  Fig.  6-4.  As  expected,  the  response  of  the 
system  becomes  faster  with  increasing  control  authority.  In  Fig.  6-5,  the  response  of  the  system  to 
a  center-span  transverse  impulse  is  shown.  This  figure  can  be  compared  with  the  open-loop  case, 
shown  in  Fig.  5-3.  In  this  case,  most  of  the  disturbance  has  been  suppressed  before  it  reflects 
from  the  boundaries,  and  the  first  mode  of  vibration  is  never  established.  This  is  characteristic  of 
high-gain  control  systems,  which  provide  high  levels  of  damping  augmentation.  For  these 
systems,  the  energy  in  the  propagating  wave  is  effectively  absorbed  as  it  progresses  towards  the 
boundaries  of  the  system, 

Li  order  to  compare  the  performance  of  the  optimal  distributed  controller  with  the 
performance  of  discrete  controllers,  finite  element  models  of  the  beam  system  were  developed. 
These  models  have  as  state  variables  the  same  quantities  that  are  used  in  the  distributed  model  (i.e., 
curvature  and  velocity),  but  are  only  available  at  discrete  points  along  the  struemre.  Similarly,  the 
control  inputs  are  available  at  a  finite  number  of  stations  along  the  beam.  These  models  were  used 
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to  develop  full  state  LQR  control  laws  for  the  discrete  systems.  Both  the  sinusoidal  initial 
condition  and  the  center-span  impact  cases  were  addressed.  Table  6-1  summarizes  the  results  of 
the  study.  In  all  cases,  the  optimal  cost  required  to  return  the  system  to  rest  approaches  the  optimal 
distributed  cost  as  the  discretization  becomes  finer.  (Details  on  deteimining  the  optimal  cost  for  the 
disBibuted  controller  can  be  found  in  Appendix  D.)  This  is  convincing  evideisce  that  the 
distributed  control  formulation  indeed  converges  upon  an  optimal  solution. 

6.3  Distributed  Control  of  an  Infinite  Beam 

In  this  section,  we  validate  the  results  for  the  optimal  control  of  fi.niie  beams  by  considering 
an  infinite  beam  system.  Most  of  the  formulation  presented  in  this  section  is  based  on  recent  work 
by  deLuis  (1989).  The  basic  idea  is  to  work  in  the  spatial  fi-equency-domain,  using  the  spatially 
transformed  dynamics  of  the  beam  system.  This  reduces  tiie  distributed  control  problem  to  a 
family  of  conventional  optimal  control  problems,  parametrized  by  the  spatial  frequency  variable. 
The  infinite  model  requires  that  the  section  properties  and  cost  weightings  be  spatially  constant,  so 
that  the  spatial  transform  is  possible.  In  his  work,  deLuis  relies  on  functional  analysis  arguments 
to  justify  the  form  of  the  control  law.  In  contrast,  tiie  approach  presented  here  is  somewhat  more 
straightforward  and  intuitively  satisfying. 

6.3.1  Spatially  Transformed  Dynanucs  and  Cost  Functional 

Because  the  beam  is  of  infinite  e.xtent,  we  can  take  the  spatial  Fourier  transform  of  (5. 1). 
The  resulting  equation  of  motion,  expressed  in  terms  of  t  and  the  spatial  frequency,  is  then 

fi(q,t)  =  A(q)5(q.t)  +  bQ(it)  (6.54) 

where  A(q)  is  obtained  from  Eq.  (5.12a): 
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1 

Weightings 

Lusaped  Paraoster  Fornulatioo 

Distributed 

Forouiation 

B 

B 

« 

elecents 

16 

eleoeots 

1 

m 

10-* 

"  ■  n 

13.63 

13.23 

13.19 

13.15 

H 

H 

10-3 

22.02 

21.72 

21.65 

21.64 

■ 

II 

D 

10-^ 

56.79 

56.67 

56.64 

56.53  1 

■ 

H 

H 

10-- 

153.43 

156.62 

156.13 

156.13  1 

1 

B 

B 

10-* 

6.50x10'* 

6.81x10'* 

6.94x10'* 

7.07x10'*  1 

iH 

H 

H 

10-’ 

2.16x10'* 

2.24x10* 

2.21x10'* 

2.24x10'-  1 

■ 

H 

H 

10-- 

6.99x10'* 

7.03x10'* 

7.05x10'* 

7.07x10'*  1 

■ 

B 

B 

10-- 

0.153 

0.160 

0.163 

Table  6-1;-  Comparison  between  optimal  costs  for  distributed  control  and  convenrional  lumped- 
parameter  control.  In  the  table,  Jj  corresponds  to  an  inirial  displacement  field  V0(x)=sin(2!ix).  and 
corresponds  to  an  initial  unit  impulsive  disturbance  applied  at  the  center-span  of  the  beam. 
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The  transformed  sate  variable  and  the  transformed  control  input  are,  in  general,  complex-valued. 
However,  because  A©  and  b  are  purely  real  for  this  system,  the  real  and  imaginary  dynamics 
decouple,  giving 

j,(5.t)  =  A©^,(5.0  +  bfl.(it) ,  ^i(tt)  =  A©<ii(q,t)  +  bO.(?.t)  (6.56a,b) 

Thus,  for  each  value  of  we  are  left  with  two  identical  real-valued  systems  for  which  an  opdmal 
control  solution  is  desired. 

The  cost  functional  to  be  minimized  is  quadratic  in  the  normalized  variables.  For 
distributed  control,  we  must  integrate  over  the  entire  (infinite)  domain.  Tne  cost  functional  is  thus 


In  this  last  expression,  which  is  analogous  to  Eq.  (6.21),  q^  and  q^  weigh  potential  and  kinetic 
energy,  respectively,  while  r  weighs  conuol  effort.  Equation  (6.57)  can  be  written 


J 


(6.5S) 


where 


^^(0  =  J[)(x,t)^y(x,t)-i-ru(x,t)i]  dx  (6.59) 

and  the  following  defutidons  have  been  made: 

y(x,t)  =  Q''^x(x,t) ,  Q  =  (6.60a,b) 

By  making  use  of  Parseval’s  theorem,  we  can  write 

•o 

J4(0  =  ^  j  [  5(5.0”Q^(^t)  +  r  Ififit)!'']  d^  (6.61) 
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where  the  superscript  represents  the  complex  conjugate  transpose  operation.  Now, 
interchanging  the  order  of  integration  in  (6.58)  yields 

eo 

J  =  i/J,(?)d5  (6.62) 

where 

00 

Ji(?)  =  j  [^(?.t)"Qx(q,t)  +  r  Iu(5,t)|^]  dt  (6.63) 

At  this  point,  the  following  observation  can  be  made:  J  is  minimized  if  and  only  if  Ji(;)  is 
minintized  for  every  value  of  We  must  now  express  and  0(^,0  in  tenns  of  their  real  and 
imaginary  parts.  Because  Q  and  r  are  purely  real,  the  (imaginary)  cross-terms  cancel  in  (6.63), 
and  we  are  left  with 

»o 

Jt(4)  =  J  [^,(it)«QV?,t)  +  r^,(;.t)2]  dt 
0 

e« 

+  I  [^i(5.0”Qx,(q,t)  +  r  0-]  dt  (6.64) 

0 


6.3.2  Optimal  Control  Solution 

For  each  value  of  q,  we  have  two  identical  dynamic  systems  given  by  Eqs.  (6.56a, b)  and 
two  identical  cost  functionals  given  by  Eq.  (6.64).  Therefore,  the  control  laws  relating  0,  to  and 
fij  to  will  be  identical,  and  can  be  combined  into  the  single  equation 

(i(q,t)  =  -(t(q)'’^^i(^,t)  (6.65) 

where  the  transformed  gain  vector,  is  a  real  function  of  the  spatial  frequency.  Taldng  the 
inverse  transform  of  this  equation  yields 
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(6.66) 


u(x.t)  =  -  Jk(x-y)'’'x(y,t)  dy 


■Dius,  the  integration  kernel  k(x-y)  can  be  though  of  as  a  weighting  from  the  sensed  state  at  a 
location  y  to  the  control  actuation  at  a  location  x.  Classical  LQR  theory  gives,  as  the  optimal 
solution 


=  [fei(§)  lc2(?)]  =  i  b'^sa)  (6.67) 

where  8(5)  solves  the  control  algebraic  Riccad  equation 

■^®^S(q)  +  S(q)A(q)  +  Q-|s(§)bb^S(q)  =  0  (6.6S) 

Substituting  the  known  parameters  A(^).  b.  Q  and  r  and  solving  for  die  elements  of  S  gives 

®u(?)  *=  l+^u’^Xf  +  2 f ■\/7+Xu  - 1 J 
Sl2(?)  =  -r§2[Vl+Xi;- 1] 

S22(?)  =  r?-'s/^T  +  2['\/T+Xi,-l] 


(6.69a) 

(6.69b) 

(6.69c) 


where 


(6.70a, b) 

Note  that  the  entire  behavior  of  the  solution  is  parametrized  by  the  dimensionless  groupings,  Xy 
and  Xj.  Since  ^  has  the  units  of  inverse  length,  we  can  infer  that  (r/qu)i/^  and  (r/qT)l«  represent 
nondimensional  distances. 

Substituting  Eq.  (6.69b)  in  Eq.  (6.67),  we  obtain,  for  the  feedback  gain  relating  curvature 
to  force. 


J=i(x)  =  (quA)^/**  fiffqu/r)'^] 
where  f,(.)  is  the  inverse  transform  of  f  j(*),  given  by 


(6.71) 
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m  =  -  + 1  - 


(6.72) 


A  plot  of  fi(*)  is  shown  in  Hg.  6-6.  Some  qaalitarive  features  of  the  feedback  gain  become 
apparent  upon  examination  of  Eq.  (6.71).  Rrst,  the  magnitude  of  the  feedback  varies  as  (qu/r)^/'*, 
so  that  increased  curvature  penalty  and  reduced  control  penality  both  increase  the  feedback  gain,  as 
expected.  Second,  the  argument  of  fi(*)  indicates  that  the  control  becomes  more  localized  with 
increasing  state  penalty  and  increasing  control  authority.  This  makes  sense  in  terms  of  the 
nondimensional  lengths  described  above.  Pfigh  control  authority  suggests  that  a  disturbance  can  be 
suppressed  quickly,  before  the  majority  of  the  energy  travels  very  far  along  the  beam,  whereas  low 
authority  requires  a  longer  time  interval  (and  hence  greater  distance)  to  suppress  the  disturbance. 
These  feamres  are  also  observed  for  the  finite  beam  system  described  in  Section  6.2.  In  fact,  cross 
sections  of  the  finite  beam  gain  kernels,  taken  near  the  center  of  the  surface,  have  the  approximate 
shape  of  the  gain  kernel  for  die  infinite  beam  system,  and  the  approximation  gets  better  with 
increasing  control  authority.  A  quantitative  analysis  indicates  that  this  is  the  case  when  r/q^  <  lO’^. 

In  computing  k2(x),  the  velocity  to  force  feedback  term,  an  interesting  feature  emerges. 

The  Riccati  solution,  S2i(\),  does  not  go  to  zero  as  ^  approaches  infinity.  As  a  result,  the  inverse 
transform  of  ^2©  include  a  delta  function.  In  order  to  make  the  inverse  transfotm 
continuous,  this  bias  term  is  subtracted  from  ic2(q),  and  a  delta  function  with  magnitude  equal  to 
this  bias  is  added  to  k2(x)  after  inversion.  Thus,  the  velocity  feedback  gain  kernel  is  expressed  as 

k2(x)  =  (qu/r)2^^  f2[(qo/>')'/‘‘x ;  qy/qu]  + 

where 

f2(5 :  Y)  =  VY-2k2fj(^)  .  VriT  (6.74) 

This  corresponds  exactly  with  the  introduction  of  a  collocated  velocity  feedback  term  for  the  finite 
beam  system.  Indeed,  the  magnitude  of  the  collocated  gain  agrees  with  the  finite  case,  with  the 
assumption  of  constant  section  properties  and  weightings.  Note  that  the  velocity  feedback  term  is 
paiametrized  by  the  same  nondimensional  length  as  ki(x),  but  an  additional  parameter,  the  ratio 
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between  kinetic  and  potential  energy  penalties,  is  also  present  Plots  of  f2(* ;  i),  shown  for  various 
y,  are  shown  in  Fig.  6-7.  Once  again,  cross  sections  of  the  velocity  feedback  gain  kernel  for  the 
finite  beam  case  agree  quite  well  with  the  infinite  beam  kernel. 

6.3.3  The  Case  of  Curvature  Actuation 

We  now  study  the  case  where  a  distributed  actuator  capable  of  inducing  curvamre  in  a 
continuous  manner  represents  the  control  input  Such  is  the  limiting  case  of  a  beam  with  many 
embedded  piezoelectric  actuators  distributed  along  its  span.  For  this  system,  the  equation  of 
motion  is  modified  to 

^v(.'(.t)  +  ^v(.x,t)  =  Ij^mefx.t)  (6.75) 

where  tri;(x,t)  represents  the  net  action  of  the  distributed  piezoelectrics.  By  ma-king  the  new 
definitions 

u(x.t)  =  m<;(x,t) ,  b(q)  =  (6.76a, b) 

the  transformed  equation  of  motion  becomes 

^(;.t)  =  A(q)^(q.t)  +  b(;)fi(;,t)  (6.77) 

This  equation  is  identical  to  Eq.  (6.54)  except  that  b  is  now  a  function  of  q.  This  subtle  difference 
has  a  profound  effect  on  the  control  law.  Following  the  same  procedure  as  in  the  previous 
subsection  leads  to  expressions  for  the  transformed  gain  kernels  which  are  independant  of  §; 

15,  =  -  1 ,  lc2  =  •\/^+2lc,  (6.78a.b) 

As  a  result,  we  have 

kiCx)  =  iczCx)  =  lcz6(x)  (6.79a,b) 
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The  optimal  control  is  therefore  purely  collocated,  even  though  the  controller  has  access  to  the  state 
vector  over  the  entire  spatial  domain.  This  result  is  quantitatively  idendcal  to  the  finite  beam  case 
presented  in  Section  6.2.5,  regardless  of  the  state  and  control  effort  penalties. 

6.4  Discussion 

It  becomes  clear  that  the  choice  of  actuator  for  the  Bemoulli-Euler  beam  system  has  a 
profound  effect  on  the  optimal  control  law.  With  a  distributed  force  acmator,  the  curvature 
feedback  is  purely  distributed  (becoming  more  localized  with  increasing  control  authority),  while 
the  velocity  feedback  has  both  collocated  and  distributed  components.  For  a  finite  beam,  a 
nondimensional  length,  which  depends  on  the  state  and  control  effort  penalties,  indicates  the  extent 
to  which  the  boundary  conditions  imposed  on  the  finite  beam  affect  the  optimal  control  solution.' 
Numerical  examples  for  finite  beam  systems  support  this  claim. 

For  a  beam  with  a  distributed  curvamre  actuator  (a  more  realistic  and  implementable 
situation),  both  the  curvature  and  velocity  feedback  gains  are  purely  collocated,  regardless  of  the 
nondimensional  length  parameter.  As  a  result,  the  boundary  conditions  do  not  affect  the  optimal 
control  solution  for  this  type  of  actuation.  The  next  logical  step  in  this  analysis  would  be  to  study 
the  effect  of  replacing  the  distributed  actuator  with  a  set  of  discrete  controls,  which  better  reflects  a 
physically  realizable  controlled  structure.  It  would  be  interesting  to  observe  whether  or  not  the 
optimal  feedback  gains  are  still  collocated  fora  set  of  discrete  embedded  piezoelectric  actuators. 

At  present,  no  claims  can  be  made  co.nceming  the  robusmess  of  the  distributed  controller. 

A  quantitative  robusmess  analysis  would  help  determine  the  sensitivity  of  the  performance  of  the 
system  to  errors  in  the  structural  model.  Also,  the  assumption  of  a  truly  distributed  controller  is 
rather  restrictive.  Any  implementable  system  will  consist  of  a  finite  set  of  sensors  and  actuators. 
Consequently,  the  theory  must  be  e.xtended  to  account  for  discrete  sensing  and  actuation.  It  may 
be  possible  to  extend  the  optimal  output  feedback  approach  discussed  by  Levine  (1971)  or  the 
optimal  projection  approach  developed  by  Bernstein  (1986)  in  a  manner  amenable  to  the  discrete 
sensing/acroation  problem. 
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Another  future  research  topic  is  the  deteimination  of  optimal  disaibuted  conaol  laws  for 
arbiaaiy  boundary  conditions.  Such  a  development  was  anempted  by  Tzafestas  (1970),  but  has 
found  exaemely  limited  application.  For  example,  the  optimal  distributed  conaol  of  a  candlevcred 
beam  was  addressed  by  Bailey  (1984)  using  Tzafestas'  formulation,  with  the  conclusion  that  the 
boundaiy  conditions  could  not  be  posed  in  the  specific  mathematical  form  required  by  the 
fonnulation.  Clearly,  the  problem  lies  in  dealing  with  the  boundaiy  conditions  which  arise  when 
determining  the  adjoint  operators  in  Eqs.  (6.6a-c).  These  boundaiy  terms  result  in  additional 
necessary  conditions  for  optimality,  expressed  in  terms  of  ordinaiy  differential  equations. 
Currently,  no  general  formulation  exists  which  includes  these  extra  conditions.  The  ability  to 
handle  general  boundaiy  conditions  would  make  it  possible  to  develop  control  laws  for  multiple 
element  structures,_sueh  as  space  frames  and  trusses. 

Another  possible  application  of  distributed  control  theory  is  in  the  active  control  of  two- 
dimensional  structures,  such  as  mirror  surfaces  and  shell  structures.  However,  numerical 
solutions  for  plates  and  membrane  models  require  extensive  computational  capabilities,  and 
therefore  represent  an  ambitious  undertaking. 
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7.  HYBRID  MODELLING  AND  CONTROL  APPROACH  FOR  HAC/LAC  DESIGN 
It  remains  to  develop  a  control  strategy  that  utilizes  the  best  aspects  of  both  the  TEM  and 
direct  modelling  methodologies.  For  example,  the  distributed  control  solutions  obtained  through 
the  direct  model  could  form  the  basis  for  a  LAC  design.  The  resulting  model  of  the  controlled 
structure  could  then  be  transformed,  and  a  HAC  controller  could  be  designed  by  posing  the 
problem  in  the  standard  form,  as  described  in  Section  4.2.  Finally,  command  prefiltering  of 
control  uiputs  for  slew  maneuvers  would  be  determined  using  the  open-loop  optimal  control  theory 
discussed  in  Section  4.1.  The  exactness  of  the  theory  makes  it  more  attractive  than  modal-based 
approaches,  such  as  the  work  of  Singer  (1990).  The  entire  hierarchically  controlled  system  would 
then  have  the  general  form  shown  in  Fig.  7-1. 

Inherent  in  this  objective  is  a  general  unification  of  the  ^vo  modelling  approaches,  which 
has  not  been  achieved  to  date.  Such  a  unification  would  be  a  profound  improvement  in  the  alilly 
to  develop  exact  control  models  for  large  flexible  structures.  Analytic  TEM  solutions  do  exist, 
however,  for  some  specific  controlled  structural  elements.  Consider,  for  example,  the  Bernoulli- 
Euler  beam  with  curvature  feedback.  The  dynamics  equation,  e.xpressed  in  dimensional  form,  is 

El|^v(x,t)  +  pAUv(x,t)  =  |^m„(x,t)  (7.1) 

The  optimal  distributed  controller  is  collocated  in  this  case,  with  the  normalized  feedback  law  given 
by 

rau(x,t)  =  -k,|^v(.x,t)  -  k2|v(x,t)  (7.2) 

where  kj  and  k2  are  determined  via  Eqs.  (6.53a,b).  Converting  this  feedback  law  into  dimensional 
form  yields 

mu(x,t)  =  -k,EI  |^v(x,t)  -  k2‘'/pAEI^v(x.t)  (7.3) 
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and  substituting  this  expression  in  the  dynamics  equation  produces  the  equation  of  motion  for  the 
controlled  structural  element.  It  is  given  by 


EI(l+ki)  ^v(x.t)  +  IC2Vp^3|ajv(x,t)  +  pA^v(x,t)  =  0  (7.4) 


It  is  now  possible  to  transform  the  closed-loop  dynamics  into  the  frequency-domain.  Assuming 
zero  mitial  conditions,  we  have 


(l+ki)|^v(x.s)  -  jk2a2|^v(x.s)  +  a^v(x,s)  =  0 


5^  :,f. 


(7.5) 


where 


=  -%r  s^ 


£I 


(7.6) 


The  homogeneous  Solution  vector  is  tlien  given  by 


The  expression  for  the  internal  state  vector  in  terms  of  v(x,s)  remains  unchanged,  and  is  given  by 
Eq.  (2.60).  The  same  can  be  said  for  the  generalized  boundary  displacements  and  forces.  The 
homogeneous  solution  vector  can  then  be  used  to  derive  analytical  expressions  for  ahe  dynamic 
stiffness  and  inteipolation  matrices,  which  will  be  slightly  more  complex  than  the  matrices 
corresponding  to  the  uncontrolled  beam. 

Unfominately,  if  the  feedback  is  indeed  distributed  rather  than  purely  collocated,  Eq.  (7  5) 
becomes  an  integral  equation.  A  general  solution  is  therefore  not  avmlable.  However,  the  form  of 
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ihe  gain  kernels  for  a  particular  problem  (e.g.,  beam  with  force  feedback)  may  lead  to  some  form 
of  anal>’ncal  solution,  or,  at  least,  an  accurate  approximate  solution.  Whether  or  not  these  classes 
of  controlled  smictures  lend  themselves  to  analytic^  TEM  models  in  the  general  case  remains  an 
unresolved  issue. 
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8..  CONCLUSIONS  AND- RECOMMENDATIONS,*!  r'nL  . ;  ■; 

The  mathemsdcally  exact  TEM  and  direct  smicmra]  modelling  methods  have  been- 
developed  and  demonstrated.  By  retaining  the  dynamics  that  describe  the  stnictura!  model  over  the 
entire  frequency  range,  accurate  behavioral  predictions  are  available.  For  example,  the  wave-like 
propagation  chaiacterisrics  associated  with  impulsive  disturbances  are  easily  observed  using  either 
exact  modelling  approach.  In  addition,  these  models  do  not  require  modal  analysis  techniques, 
although  modal  information  is  available  for  TEM  models  of  frame-like  structures.  Furthermore, 
the  frequency-domain  analysis  incorporates  general  viscoelastic  damping  mechanisms  in  a  very 
straghtforward  manner.  Rnally,  the  dramatic  increase  in  computation  speed  of  the  TEM  artalysis 
technique  over  traditional  finite  element  modelling  has  been  demonsuated. 

New  control  formulations  were  developed  that  take  advantage  of  the  infomation  available 
via  the  exact  modelling  methods.  An  open-loop  optimal  control  technique  was  denaonsaated  using 
TEM  models,  and  was  found  to  virtually  eliminate  the  residual  enerjy  associated  with  the  slewing 
of  flexible  struentres.  The  only  approximation  made  concerned  the  control  inputs  themselves, 
w'lu'ch  were  assumed  to  have  lirtuted  bandwidth.  The  direct  analysis  technique  provided  the 
framework  for  a  distributed  control  theory.  Having  been  developed,  the  theory  was  applied  to  a 
simple  Bemoulli-Euler  beam  system,  and  the  feedback  gain  kernels  were  determined.  These 
kernels  agreed  with  previous  results  concerning  the  optimal  distributed  control  of  an  inEiite  beam. 

Several  issues  remain  unresolved,  and  are  recommended  for  future  research.  Although  the 
TEM  methodology  has  been  developed  for  two-dimensional  structures,  it  is  incomplete  in  two 
respects.  First,  the  selection  of  boundary  points  and  their  relation  to  element  geometry  and 
solution  accuracy  must  be  addressed.  Second,  a  rigorous  method  of  determining  the  set  of  basis 
solutions  for  the  homogeneous  solution  vector  must  be  developed.  Clearly,  these  two  goals  are 
intimately  tied,  as  the  number  of  basis  solutions  required  is  directly  related  to  the  number  of  the 
boundary  points  used.  Inherent  in  this  analysis  is  a  comparison  between  the  TEM  analysis  and 
finite  element  models  of  two-dimensional  structures. 
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^  x'rM  The  cfosed-Ioo?  control  problea,  expressed  in  the  fiequency-domain  with  TEM  models, 
remains  to  be'  solved?-  Here;  the  lack  of  a  finite  state-space  represenation  of  the  plant  is  the- 
fundamental  difficulty.  It  may  be  possible  to  extend  the  coprime  factorization  technique  to  general 
strucwral  systems,  or  develop  some  other  approach. 

The  distributed  control  solutions  developed  here  must  be  extended  to  include  other 
structural  elements,  such  as  Timoshenko  beams  and  axial  and  torsional  rods.  Two-dimensional 
elements,  which  may  represent  deformable  mirror  surfaces  or  solar  panels,  must  also  be 
incoiporated  into  the  distributed  control  framework,  although  this  extension  presents  a 
considerably  more  difficult  computational  challenge.  Finally,  the  theory  roust  be  extended  to  apply 
to  multi-element  structures,  which  may  include  any  of  the  elements  mentioned  The  ability  to 
handle  the  complex  boundary  conditions  that  arise  at  element  junctions  is  tiie  primary  difficulty 
here. 

An  evaluation  of  the  robustness  of  the  distributed  controller  to  model  uncenainty  must  also 
be  undertalten.  Discrepencies  between  the  model  and  the  actual  physical  struemre,  caused  by 
tolerances  in  physical  dimensions  and  material  properties,  structural  joint  dynamics,  nonlinear 
material  behavior  and  other  unmodelled  dynamics,  usually  result  in  perforaiance  degradation.  A 
rigorous  robusmess  evaluation  would  quantify  the  relation  between  modelling  error  and 
performance.  Linked  to  this  issue  are  implementation  considerations.  Because  actuators  and 
sensors  are  always  discrete  in  nature,  the  distributed  control  solution  represents  only  a  limiting  case 
as  the  number  of  individual  actuators  and  sensors  approaches  infinity.  For  any  '.nplementable 
control  design,  then,  the  effect  of  utilizing  a  finite  set  of  actuators  and  sensors  must  be  addressed. 
Also  relevant  is  a  study  of  the  effect  of  actuator  and  sensor  dynamics  and  their  relation  to  robust 
stability. 

A  combined  direct/TEM  control  methodology  for  structural  systems  is  not  yet  available. 
This  hybrid  technique  would  facilitate  the  development  of  hierarchical  control  schemes,  such  as 
HA  C/LAC,  without  resorting  to  modal  analysis  and  truncation.  Consequently,  the  problems  of 
control  and  observation  spillover  would  be  alleviated,  at  least  firom  a  mathematical  perspective. 


conservative  control  design  would  be  required,  resulting  in  enhanced  notnimal  performance.- 
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APPENDK^Ar  HIGH  FREQUENCY,  TEM  ELEMENTS 

"  '^1*  « 

This  ap^hdix  presents  the  dyniiiic  stiffnes  matrices  for  high  frequency  TEM  elements. 
In  pardculaf,  *e  Mindlin-Hemnann  axial  rod  and  she  Timoshe^o  be^,  Ixith  described  in  Sec. 
2.3.6,  are  discussed. 


A.l  Mindlin-Herrmann  Rod 

The  dynamics  of  the  Mindlin-Herrmann  rod  are  characterized  by  the  following  set  of 
differennal  equanons: 

-a2(X+2G)|^vi(x,s)  '•paVvifx.s)  =  2aX|^V2(x,s)  (A.la) 

a2K2G^V2(x.s)-  {^8K^(X+G)  +  pa2s2jv2(x,s)  =  4aK^x|j  vi(x,s)  (A. lb) 

All  symbols  in  these  equadons  are  defined  in  Sec.  2.3.6.  These  equations  apply  to  a  rod  of 
circular  cross-secrion  only.  The  Lame  constants  are  related  to  the  modulus  of  elasticity  and 
Poisson's  ratio  of  the  material  according  to 


(!+v)(l.2v) 


.  G  = 


(A.2a,b) 


VVe  will  assume  that  the  radial  deformation,  V2,  is  constrained  at  the  boundaries  of  the  element 
This  is  the  case  when  the  element  is  embedded  in  a  stnictuial  junction.  We  therefore  have 

V2(0,s)  =  V2(L,s)  =  0  (A 

With  these  constraints  imposed,  the  stiffness  matrix  becomes  2-by-2,  and  can  be  expressed  by 


KM  -  •PiSt  +  p2S2  1 

^  ’  ~  ri(s)  L  -  PiSj  +  P2S2  PlSi'2  ■  P2S2':i J 


am  -  ('^PlP2(>-Cl<=2)  (Pi^P^SiSa] 

(1-v)  L  Pl(“2a)-P2(Oia)  J 


and  the  following  trigonometric  definitions  have  been  made: 
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VT,  HOin  :A 


'i,*""  ^lii  'I'cy* 


where 


1  _  tlii£!^  •»  _  n+x^UI-2v^  4-  1 

^1  “  v2/i.>,\  ^2  “  „i 


X3  = 


x2(l.v) 

_  (UvVl-2v) 


Sk; 


pa£s2. 


4x: 


(A.6a,b) 


=-^V^"o~' J a ’V  5--“--  trr^I-xr c-.i.tr?.r!-.'^L':';/:s.'i'.--rrt-,r  'j.-J 

Also,  Pi  and  P2  are  given  by 

^  [ic2(l+v)(l-2v)2pa2s2.i6ic^2El((x.a)4.x2(l.v)(J-2v)E(a,a)2 

Pi  —  - — ir2 - ?  (A./) 

‘  4(l-v)(l-2v)paV  +  l6xJE 

and  the  nondimensional  paiaineters  aia  and  034  are  given  by 

(aia)2  =  -X,[(l+X2a)±-\^(1+X2a)2-5^(1+X30)a  (A.8) 


(A.9a-d) 


The  stiffness  matrix  can  be  shown  to  reduce  to  that  of  the  simple  axial  rod  either  by  setting  v=0  or 
by  taking  the  limit  as  the  radius,  a,  approaches  zero. 


A.2  Timoshenko  Beam 

The  dynamics  of  the  Timoshenko  beam  can  be  expressed  either  as  a  system  of  two  coupled 
partial  differential  equations  or  as  the  single  equation 

EI^v(x,t)  +  pAv(x,t)  -  pi  [l-t~]|^v(x,t)  +  ^v(x,t)  =  f(i(x,t)  (A.IO) 

Taking  the  Laplace  transform  of  this  equation  yields 

^v(x,s)  +  2pia20v(x,s)  -  a^v(x,s)  =  f(i(x,s)  (A.ll) 

where  the  following  paranwters  have  been  defmed: 
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(A.12a-e) 


Pi.=  5(l+ki)(ar)? .  p|  =  i -ki(ar)4,. 


I 

A 


The  homogeneous  solution  vector  is  then 


(50  -  f.’' 


Vh(x,s) 


■cosh(apix)' 
sinh{apix) 
cos(aP2x) 
.  sin(aP2x)- 


where 

P?  =  P3-Pl.  P2  =  P3  +  P1 
and 

P3  =  P1+P2  =  l+|(l-ki)-(ctr)^ 


(A.13) 


(A.14a.b) 


(A.I5) 


Due  the  the  internal  shearing  allowed  by  the  'Ilnioshenko  model,  the  expressions  for  the  internal 
state  vector  in  terms  of  the  homogeneous  solution  become  considerably  more  complex.  They  are: 


where 


u(x,s) 


rv(x.s)l 

6(x,s) 

M(x.s) 

Ls(x.s)J 


(A.16) 


P4  =  ici(ctr)2 


(A.17) 


The  stiffness  matrix  then  takes  the  form  given  by  Eq.  (2.66),  with 
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ACTUATION 


In  this  appendix,  we  extend  die  results  of  Secdon  5. 1.2.3  to  the  case  of  a  beam  with  a 
distributed  curvature  actuator.  The  modified  equation  of  motion,  expressed  in  the  time-domian,  is 

mu(x,t)  +  f„(x,t)  (B.l) 

which  becomes,  after  Laplace  transformation 

0  riiuCx.s)  +  f„(x,s)  (B.2) 

The  assumed  linear  feedback  control  law  is 
1 

mu(x,0  =  -j  [ki(x.y)Ti(y)^v(y,i)  +k2(x,y)|v(y,t)j  dy  (B.3) 

which  transforms  into 


mu(x.s)  =  •J[l<i(x,y)T|(y)^v(y,s) +k2(x,y)sv(y.s)jdy 

1 

+ J  [k2(x,y)vo(y)]  dy 


(B.4) 


Substituting  this  expression  into  Eq.  (B.2)  leads  to 

1 


+  skaCx.y)  v(y,s)j  dy 


fn(x.s)  +  ^  [vo(x)  +  s  vo(x)]  +  ^ 


ii 


lc2(x,y)  vo(y)  dy 


(B.5) 


Employing  the  same  discretization  technique  used  in  Section  5. 1.2,3,  we  obtain,  for  the  finite- 
difference  equation 


1 

1 

i 


lae 


2  J. 


LSI n^KX?)]  V§X^  NjD ?<;($).,,  .  ^  (B.6) 


All  terms  in  this  last  equation  are  as  defined  in  Section  5.1.Z3. 


147 


'*  appendix;  Cj!  SmUtATION  OP  A.TmOSHENKO  BEAM* 

-  *■  -  In‘Wsapp«ndix.wediscussthe(ijrectsimtiIanonofarm»sbenJcobemF6rM^^^^^^ 

we  will  assumed  that  the  initial  conditions  are  zero  and  that  the  secdon  piopetties  are  constant, 
although  other  situarions  can  be  treated  in  a  manner  sinular  to  the  Bemoulli-Euler  beam  model. 

The  dimensional  form  of  the  equation  of  motion  is 

EI$v<,(Xd.t<,)  +  pA^v,(Xd.t,)  -  pi  [i4f]^v„(x<j.ti)  +  f=^v,(x,,t,) 

=  f<i(>:.0  (C.1) 

which  can  be  normalized  (using  the  same  groupings  as  in  Section  5.1.2.1)  to 

^v(x,t)-+  |jv(x.t)  -  a,[i+a£]^v(x,t)  +  ajJag v(x.t)  =  f(x,t)  (C2) 


~  I  Ek 

=  «E  =  f 


(C.3a.b) 


Talang  the  Laplace  transform  (neglecting  initial  conditions)  jields 

§v(x.s)  -  o:,[l+ctE]s2|^v(x.s)  +  s2v(x,s)  +  ai2aEs^v(x,s)  =  f(x,s)  (C.4) 
Finally,  the  spatial  discretization  yields 

[N^  D 2  +  s2  (I  -  0t,(l+aE)N'2 D)  +  a,2aE  s^lj  v(s)  =  f(s)  (C.g) 


/  ;<• 
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APPENDIX  DV  O^tlMAL'’ COSTS  FOR  DISTRIBUTED  CONTROL’ SYSTEMS 


1  This  appendix  discusses  die  computation  of  the  optimal  cost  reqmied  to  bring  a  distributed 
system  to  rest  from  an  arbitrary  inidal  condition  using  a  distributed  controller.  We  will  restrict 
attention  to  the  force  actuation  case,  and  we  will  assume  uniform  cross  section  properties.  The 
formulation  is  analogous  to  the  discrete  case,  for  which  the  optimal  cost  is  expressed  by 

r  =  |xJSxt,  (D.l) 

where  Xg  is  the  initial  condition  on  the  state  vector  and  S  is  the  Riccati  matrix.  Wang  (1964) 
shows  that,  for  the  distributed  case,  the  optimal  cost  has  the  foim 

1 1 

J*  =  iJJxo(x)'''S(x,y)Xo(y)dxdy  (D.z) 


where  Xo(x)  is  the  distributed  initial  condidon  and  S(x,y)  is  the  Riccad  matrix  function  associated 
with  Eq.  (6 19).  Thus,  for  the  case  of  an  inidal  displacement,  vo(x)=sin(2mc),  the  optimal  cost  is 


J*=i 


11,  j 

2  JJ  ^''oW  Sii(x.y)  |;2''o(y)  dx  dy 


1 1 


=  87c4  /Jsin(27tx)  Sji(x,y)  sin(2xy)dxdy 


(D.3) 


The  Riccad  equations  for  Sii(x,y)  follow  direcdy  from  Eq.  (6.25).  For  constant  section 
properties,  they  are  given  by 
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3;;2^u(*>y)  '  .  i J S22(x,z)Si2(z,y)dz  ^  O’' 

'V'  >  o  J'.,  i  '  r  .  '  t  '  '  * 

32  32  ^ 

*  a^2S22(x.y)  +  3;2Sii(x.y)  -  ^  JSi2(x.z)Si2(z.y)dz  =  0  J  (D.4b) 
Adding  these  two  equations  yields 

V2S„(x,y)  =  V2s,j(x.y)  +  j  J  [S22(x.z)S,2(z,y)^S,2(x,z)S22(z.y)Jdz  (D.5) 

Thus.  S,i  is  computed  using  a  simple  relaxation  algorithm.  For  the  case  of  an  initial  velocity, 
vo(x)=S(x-3/2),  which  is  equivalent  to  a  center-span  impulse,  the  cost  is  simply 

1 1 

“  2oy''o(’‘)S22(x,y)vo(y)dxdy  =  1/2)  (D.6) 
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